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Preface 



The purpose of this book is to present detailed fundamental information on a 
global positioning system (GPS) receiver. Although GPS receivers are popu- 
larly used in every-day life, their operation principles cannot be easily found 
in one book. Most other types of receivers process the input signals to obtain 
the necessary information easily, such as in amplitude modulation (AM) and 
frequency modulation (FM) radios. In a GPS receiver the signal is processed 
to obtain the required information, which in turn is used to calculate the user 
position. Therefore, at least two areas of discipline, receiver technology and 
navigation scheme, are employed in a GPS receiver. This book covers both 
areas. 

In the case of GPS signals, there are two sets of information: the civilian 
code, referred to as the coarse/acquisition (C/A) code, and the classified mil- 
itary code, referred to as the P(Y) code. This book concentrates only on the 
civilian C/A code. This is the information used by commercial GPS receivers 
to obtain the user position. 

The material in this book is presented from the software receiver viewpoint 
for two reasons. First, it is likely that narrow band receivers, such as the GPS 
receiver, will be implemented in software in the future. Second, a software 
receiver approach may explain the operation better. A few key computer pro- 
grams can be used to further illustrate some points. 

This book is written for engineers and scientists who intend to study and 
understand the detailed operation principles of GPS receivers. The book is at 
the senior or graduate school level. A few computer programs written in Matlab 
are listed at the end of several chapters to help the reader understand some of 
the ideas presented. 

I would like to acknowledge the following persons. My sincere appreciation 
to three engineers: Dr. D. M. Akos from Stanford University, M. Stockmaster 
from Rockwell Collins, and J. Schamus from Veridian. They worked with me 
at the Air Force Research Laboratory, Wright Patterson Air Force Base on the 
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design of a software GPS receiver. This work made this book possible. Dr. 
Akos also reviewed my manuscripts. I used information from several courses 
on GPS receivers given at the Air Force Institute of Technology by Lt. Col. 
B. Riggins, Ph.D. and Capt. J. Requet, Ph.D. Valuable discussion with Drs. 
F. VanGraas and M. Braasch from Ohio University helped me as well. I am 
constantly discussing GPS subjects with my coworkers, D. M. Lin and V. D. 
Chakravarthy. 

The management in the Sensor Division of the Air Force Research Labo- 
ratory provided excellent guidance and support in GPS receiver research. Spe- 
cial thanks are extended to Dr. P. S. Hadom, E. R. Martinsek, A. W. White, 
and N. A. Pequignot. I would also like to thank my colleagues, R. L. Davis, 
S. M. Rodrigue, K. M. Graves, J. R. McCall, J. A. Tenbarge, Dr. S. W. Schnei- 
der, J. N. Hedge Jr., J. Caschera, J. Mudd, J. P. Stephens, Capt. R. S. Parks, 
P. G. Howe, D. L. Howell, Dr. L. L. Liou, D. R. Meeks, and D. Jones, for their 
consultation and assistance. 

Last, but not least, I would like to thank my wife, Susan, for her encourage- 
ment and understanding. 




Notations and Constants 



a e - 6378137 ± 2 m is the semi-major axis of the earth. 
cifQ is the satellite clock correction parameter. 

Of i is the satellite clock correction parameter. 
a /2 is the satellite clock correction parameter. 
a s is the semi-major axis of the satellite orbit 
A bi is the satellite clock error. 

b e - 6356752.3142 m is the semi-minor axis of the earth. 
b s is the semi-minor axis of the satellite orbit 

b u is the user clock bias error, expressed in distance, which is related to the 
quantity b ut by b u = cb ut . 
b ul is the user clock error. 

c = 2.99792458 x 10 8 meter/sec is the speed of light. 

C ic is the amplitude of the cosine harmonic correction term to the angle of 
inclination. 

Ci S is the amplitude of the sine harmonic correction term to the angle of 
inclination. 

C rc is the amplitude of the cosine harmonic correction term to the orbit 
radius. 

C rs is the amplitude of the sine harmonic correction term to the orbit radius. 
c s is the distance from the center of the ellipse to a focus. 

C uc is the amplitude of the cosine harmonic correction term to the argument 
of latitude. 

C us is the amplitude of the sine harmonic correction term to the argument 
of latitude. 

A Dj is the satellite position error effect on range. 

E is called eccentric anomaly. 
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xvi 



e e - 0.0818191908426 is the eccentricity of the earth. 
e p - 0.00335281066474 is the ellipticity of the earth. 
e s is the eccentricity of the satellite orbit. 

F = -4.442807633 x 10 10 sec/m 1 / 2 . 
fj is the input frequency. 

/o is the output frequency in baseband. 
f s is the sampling frequency. 
h is altitude. 

i is the inclination angle at reference time, 
idot is the rate of inclination angle. 

A/, is the ionospheric delay error. 

I is longitude. 

L is geodetic latitude often used in maps. 

L c is geocentric latitude. 

M is mean anomaly. 

M 0 is the mean anomaly at reference time. 

An is the mean motion difference from computed value. 
r e = 6368 km is average earth radius. 

ro is the distance from the center of the earth to the point on the surface of 
the earth under the user position, 
ro, is the average radius of an ideal spherical earth. 
r s is the average radius of the satellite orbit. 
t is the GPS time at time of transmission corrected for transit time. 
t c is the coarse GPS system time at time of transmission corrected for transit 
time. 

Tgd is the satellite group delay differential. 

AT i is the tropospheric delay error. 
t oc is the satellite clock correction parameter. 
t oe is the reference time ephemeris. 
t p is the time when the satellite passes the perigee. 
t S i is referred to as the true time of transmission from satellite i. 
t t is the transit time (time for the signal from the satellite to travel to the 
receiver). 

t u is the time of reception. 
v s is the speed of the satellite. 

/x = 3.986005 x 10 14 meters 3 /sec 2 is the earth’s universal gravitational 
parameter. 

v i is the receiver measurement noise error. 

A Vi is the relativistic time correction. 




NOTATIONS AND CONSTANTS 



xvii 



7T = 3.1415926535898. 

p, 7 is the true value of pseudorange from user to satellite i. 

Pi is the measured pseudorange from user to satellite i 
o) is the argument of the perigee. 

Q f (Q - a ) is the modified right ascension angle. 

Q e is the longitude of ascending node of orbit plane at weekly epoch. 

Q er is the angle between the ascending node and the Greenwich meridian. 
U is the rate of the right ascension. 

(l ie - 7.2921151467 x 10 5 rad/sec is the WGS-84 value of the earth’s rota- 
tion rate. 
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CHAPTER. ONE 



Introduction 



1.1 INTRODUCTION) 1-13 ) 

This book presents detailed information in a compact form about the global 
positioning system (GPS) coarse/acquisition (C/A) code receiver. Using the 
C/A code to find the user location is referred to as the standard position service 
(SPS). Most of the information can be found in references 1 through 13. How- 
ever, there is much more information in the references than the basics required 
to understand a GPS receiver. Therefore, one must study the proper subjects 
and put them together. This is a tedious and cumbersome task. This book does 
this job for the reader. 

This book not only introduces the information available from the references, 
it emphasizes its applications. Software programs are provided to help under- 
stand some of the concepts. These programs are also useful in designing GPS 
receivers. In addition, various techniques to perform acquisition and tracking 
on the GPS signals are included. 

This book concentrates only on the very basic concepts of the C/A code 
GPS receiver. Any subject not directly related to the basic receiver (even if 
it is of general interest, i.e., differential GPS receiver and GPS receiver with 
carrier-aided tracking capacity) will not be included in this book. These other 
subjects can be found in reference 1. 



1.2 HISTORY OF GPS DEVELOPMENT) 1 ' 5 ' 12 ) 

The discovery of navigation seems to have occurred early in human history. 
According to Chinese storytelling, the compass was discovered and used in wars 
during foggy weather before recorded history. There have been many different 
navigation techniques to support ocean and air transportation. Satellite-based 
navigation started in the early 1970s. Three satellite systems were explored 
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before the GPS programs: the U.S. Navy Navigation Satellite System (also 
referred to as the Transit), the U.S. Navy’s Timation (TIMe navigATION), and 
U.S. Air Force project 621B. The Transit project used a continuous wave (cw) 
signal. The closest approach of the satellite can be found by measuring the 
maximum rate of Doppler shift. The Timation program used an atomic clock 
that improves the prediction of satellite orbits and reduces the ground control 
update rate. The Air Force 621B project used the pseudorandom noise (PRN) 
signal to modulate the carrier frequency. 

The GPS program was approved in December 1973. The first satellite was 
launched in 1978. In August 1993, GPS had 24 satellites in orbit and in Decem- 
ber of the same year the initial operational capability was established. In February 
1994, the Federal Aviation Agency (FAA) declared GPS ready for aviation use. 



1.3 A BASIC GPS RECEIVER 

The basic GPS receiver discussed in this book is shown in Figure 1.1. The sig- 
nals transmitted from the GPS satellites are received from the antenna. Through 
the radio frequency (RF) chain the input signal is amplified to a proper ampli- 
tude and the frequency is converted to a desired output frequency. An analog- 
to-digital converter (ADC) is used to digitize the output signal. The antenna, 
RF chain, and ADC are the hardware used in the receiver. 

After the signal is digitized, software is used to process it, and that is why this 
book has taken a software approach. Acquisition means to find the signal of a 
certain satellite. The tracking program is used to find the phase transition of the 
navigation data. In a conventional receiver, the acquisition and tracking are per- 
formed by hardware. From the navigation data phase transition the subframes 
and navigation data can be obtained. Ephemeris data and pseudoranges can be 




User Satellite lEphemeris &l Subframe 

| position I positions lpseudorangc identify 



Tracking 



J 1 Acquis ition 



Software 



FIGURE 1.1 A fundamental GPS receiver. 
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obtained from the navigation data. The ephemeris data are used to obtain the 
satellite positions. Finally, the user position can be calculated for the satellite 
positions and the pseudoranges. Both the hardware used to collect digitized data 
and the software used to find the user position will be discussed in this book. 



1.4 APPROACHES OF PRESENTATION 

There are two possible approaches to writing this book. One is a straightforward 
way to follow the signal flow shown in Figure 1.1. In this approach the book 
would start with the signal structure of the GPS system and the methods to pro- 
cess the signal to obtain the necessary the information. This information would 
be used to calculate the positions of the satellites and the pseudoranges. By 
using the positions of the satellites and the pseudoranges the user position can 
be found. In this approach, the flow of discussion would be smooth, from one 
subject to another. However, the disadvantage of this approach is that readers 
might not have a clear idea why these steps are needed. They could understand 
the concept of the GPS operation only after reading the entire book. 

The other approach is to start with the basic concept of the GPS from a 
system designers’ point of view. This approach would start with the basic con- 
cept of finding the user position from the satellite positions. The description 
of the satellite constellation would be presented. The detailed information of 
the satellite orbit is contained in the GPS data. In order to obtain these data, 
the GPS signal must be tracked. The C/A code of the GPS signals would then 
be presented. Each satellite has an unique C/A code. A receiver can perform 
acquisition on the C/A code to find the signal. Once the C/A code of a certain 
satellite is found, the signal can be tracked. The tracking program can produce 
the navigation data. From these data, the position of the satellite can be found. 
The relative pseudorange can be obtained by comparing the time a certain data 
point arrived at the receiver. The user position can be calculated from the satel- 
lite positions and pseudoranges of several satellites. 

This book takes this second approach to present the material because it 
should give a clearer idea of the GPS function from the very beginning. The 
final chapter describes the overall functions of the GPS receiver and can be 
considered as taking the first approach for digitizing the signal, performing 
acquisition and tracking, extracting the navigation data, and calculating the user 
position. 



1.5 SOFTWARE APPROACH 

This book uses the concept of software radio to present the subject. The soft- 
ware radio idea is to use an analog-to-digital converter (ADC) to change the 
input signal into digital data at the earliest possible stage in the receiver. In 
other words, the input signal is digitized as close to the antenna as possible. 
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Once the signal is digitized, digital signal processing will be used to obtain 
the necessary information. The primary goal of the software radio is minimum 
hardware use in a radio. Conceptually, one can tune the radio through software 
or even change the function of the radio such as from amplitude modulation 
(AM) to frequency modulation (FM) by changing the software; therefore great 
flexibility can be achieved. 

The main purpose of using the software radio concept to present this subject 
is to illustrate the idea of signal acquisition and tracking. Although using hard- 
ware to perform signal acquisition and tracking can also describe GPS receiver 
function, it appears that using software may provide a clearer idea of the sig- 
nal acquisition and tracking. In addition, a software approach should provide a 
better understanding of the receiver function because some of the calculations 
can be illustrated with programs. Once the software concept is well understood, 
the readers should be able to introduce new solutions to problems such as var- 
ious acquisition and tracking methods to improve efficiency and performance. 
At the time (December 1997) this chapter was being written, a software GPS 
receiver using a 200 MHz personal computer (PC) could not track one satellite 
in real time. When this chapter was revised in December 1998, the software 
had been modified to track two satellites in real time with a new PC operat- 
ing at 400 MHz. Although it is still impossible to implement a software GPS 
receiver operating in real time, with the improvement in PC operating speed and 
software modification it is likely that by the time this book is published a soft- 
ware GPS receiver will be a reality. Of course, using a digital signal processing 
(DSP) chip is another viable way to build the receiver. 

Only the fundamentals of a GPS receiver are presented in this book. In order 
to improve the performance of a receiver, fine tuning of some of the operations 
might be necessary. Once readers understand the basic operation principles of 
the receiver, they can make the necessary improvement. 



1.6 POTENTIAL ADVANTAGES OF THE SOFTWARE APPROACH 

An important aspect of using the software approach to build a GPS receiver 
is that the approach can drastically deviate from the conventional hardware 
approach. For example, the user may take a snapshot of data and process them 
to find the location rather than continuously tracking the signal. Theoretically, 
30 seconds of data are enough to find the user location. This is especially use- 
ful when data cannot be collected in a continuous manner. Since the software 
approach is in the infant stage, one can explore many potential methods. 

The software approach is very flexible. It can process data collected from 
various types of hardware. For example, one system may collect complex data 
referred to as the inphase and quadrature-phase (I and Q) channels. Another 
system may collect real data from one channel. The data can easily be changed 
from one form to another. One can also generate programs to process complex 
signals from programs processing real signals or vice versa with some simple 
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modifications. A program can be used to process signals digitized with various 
sampling frequencies. Therefore, a software approach can almost be considered 
as hardware independent. 

New algorithms can easily be developed without changing the design of the 
hardware. This is especially useful for studying some new problems. For exam- 
ple, in order to study the antijamming problem one can collect a set of digitized 
signals with jamming signals present and use different algorithms to analyze it. 



1.7 ORGANIZATION OF THE BOOK 

This book contains nine chapters. Chapter 2 introduces the user position require- 
ments, which lead to the GPS parameters. Also included in Chapter 2 is the basic 
concept of how to find the user position if the satellite positions are known. Chap- 
ter 3 discusses the satellite constellation and its impact on the GPS signals, which 
in turn affects the design of the GPS receiver. Chapter 4 discusses the earth-cen- 
tered, earth-fixed system. Using this coordinate system, the user position can be 
calculated to match the position on every-day maps. The GPS signal structure is 
discussed in detail in Chapter 5. Chapter 6 discusses the hardware to collect data, 
which is equivalent to the front end of a conventional GPS receiver. Changing the 
format of data is also presented. Chapter 7 presents several acquisition methods. 
Some of them can be used in hardware design and others are suitable for soft- 
ware applications. Chapter 8 discusses two tracking methods. One uses the con- 
ventional phase-locked loop approach and the other one is more suitable for the 
software radio approach. The final chapter is a summary of the previous chapters. 
It takes all the information in the first eight chapters and presents in it an order 
following the signal flow in a GPS receiver. 

Computer programs written in Matlab are listed at the end of several chap- 
ters. Some of the programs are used only to illustrate ideas. Others can be used 
in the receiver design. In the final chapter all of the programs related to design- 
ing a receiver will listed. These programs are by no means optimized and they 
are used only for demonstration purposes. 



REFERENCES 

1. Parkinson, B. W., Spilker, J. J. Jr., Global Positioning System: Theory and Appli- 
cations, vols. 1 and 2, American Institute of Aeronautics and Astronautics, 370 
L’Enfant Promenade, SW, Washington, DC, 1996. 

2. “System specification for the navstar global positioning system,” SS-GPS-300B 
code ident 07868, March 3, 1980. 

3. Spilker, J. J., “GPS signal structure and performance characteristics,” Navigation, 
Institute of Navigation, vol. 25, no. 2, pp. 121-146, Summer 1978. 

4. Milliken, R. J., Zoller, C. J., “Principle of operation of NAVSTAR and system char- 
acteristics,” Advisory Group for Aerospace Research and Development (AGARD) 




6 INTRODUCTION 



Ag-245, pp. 4-1—4. 1 2, July 1979. 

5. Misra, P. N., “Integrated use of GPS and GLONASS in civil aviation,” Lincoln Lab- 
oratory Journal, Massachusetts Institute of Technology, vol. 6, no. 2, pp. 231-247, 
Summer/Fall, 1993. 

6. “Reference data for radio engineers,” 5th ed., Howard W. Sams & Co. (subsidiary 
of ITT), Indianapolis, 1972. 

7. Bate, R. R., Mueller, D. D., White, J. E., Fundamentals of Astrodynamics, pp. 
182-188, Dover Publications, New York, 1971. 

8. Wells, D. E., Beck, N., Delikaraoglou, D., Kleusbery, A., Krakiwsky, E. J., 
Lachapelle, G., Langley, R. B., Nakiboglu, M., Schwarz, K. P., Tranquilla, J. M., 
Vanicek, P, Guide to GPS Positioning, Canadian GPS Associates, Frederiction, 
N.B., Canada, 1987. 

9. “Department of Defense world geodetic system, 1984 (WGS-84), its definition and 
relationships with local geodetic systems,” DMA-TR-8350.2, Defense Mapping 
Agency, September 1987. 

10. “Global Positioning System Standard Positioning Service Signal Specification, 2nd 
ed., GPS Joint Program Office, June 1995. 

11. Bate, R. R., Mueller, D. D., White, J. E., Fundamentals of Astrodynamics, Dover 
Publications, New York, 1971. 

12. Riggins, B., “Satellite navigation using the global positioning system,” manuscript 
used in Air Force Institute of Technology, Dayton OH, 1996. 

13. Kaplan, E. D., ed., Understanding GPS Principles and Applications, Artech House, 
Norwood, MA, 1996. 




Fundamentals of Global Positioning System Receivers: A Software Approach 
James Bao-Yen Tsui 
Copyright © 2000 John Wiley & Sons, Inc. 
Print ISBN 0-471-38154-3 Electronic ISBN 0-471-20054-9 



CHAPTER. TWO 



Basic GPS Concept 



2.1 INTRODUCTION 

This chapter will introduce the basic concept of how a GPS receiver determines 
its position. In order to better understand the concept, GPS performance require- 
ments will be discussed first. These requirements determine the arrangement of 
the satellite constellation. From the satellite constellation, the user position can 
be solved. However, the equations required for solving the user position turn 
out to be nonlinear simultaneous equations, which are difficult to solve directly. 
In addition, some practical considerations (i.e., the inaccuracy of the user clock) 
will be included in these equations. These equations are solved through a lin- 
earization and iteration method. The solution is in a Cartesian coordinate system 
and the result will be converted into a spherical coordinate system. However, 
the earth is not a perfect sphere; therefore, once the user position is found, the 
shape of the earth must be taken into consideration. The user position is then 
translated into the earth-based coordinate system. Finally, the selection of satel- 
lites to obtain better user position accuracy and the dilution of precision will 
be discussed. 



2.2 GPS PERFORMANCE REQUIREMENTS! 1 ) 

Some of the performance requirements are listed below: 

1. The user position root mean square (rms) error should be 10-30 m. 

2. It should be applicable to real-time navigation for all users including the 
high-dynamics user, such as in high-speed aircraft with flexible maneu- 
verability. 

3. It should have worldwide coverage. Thus, in order to cover the polar 
regions the satellites must be in inclined orbits. 
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4. The transmitted signals should tolerate, to some degree, intentional 
and unintentional interference. For example, the harmonics from some 
narrow-band signals should not disturb its operation. Intentional jamming 
of GPS signals is a serious concern for military applications. 

5. It cannot require that every GPS receiver utilize a highly accurate clock 
such as those based on atomic standards. 

6. When the receiver is first turned on, it should take minutes rather than 
hours to find the user position. 

7. The size of the receiving antenna should be small. The signal attenuation 
through space should be kept reasonably small. 

These requirements combining with the availability of the frequency band 
allocation determines the carrier frequency of the GPS to be in the L band (1-2 
GHz) of the microwave range. 



2.3 BASIC GPS CONCEPT 

The position of a certain point in space can be found from distances measured 
from this point to some known positions in space. Let us use some examples to 
illustrate this point. In Figure 2.1, the user position is on the x-axis; this is a one- 
dimensional case. If the satellite position Si and the distance to the satellite x\ 
are both known, the user position can be at two places, either to the left or right 
of Si. In order to determine the user position, the distance to another satellite 
with known position must be measured. In this figure, the positions of S 2 and 
X 2 uniquely determine the user position U. 

Figure 2.2 shows a two-dimensional case. In order to determine the user 
position, three satellites and three distances are required. The trace of a point 
with constant distance to a fixed point is a circle in the two-dimensional case. 
Two satellites and two distances give two possible solutions because two circles 
intersect at two points. A third circle is needed to uniquely determine the user 
position. 

For similar reasons one might decide that in a three-dimensional case four 
satellites and four distances are needed. The equal-distance trace to a fixed point 
is a sphere in a three-dimensional case. Two spheres intersect to make a circle. 
This circle intersects another sphere to produce two points. In order to determine 
which point is the user position, one more satellite is needed. 



£ 



s, 



t 



1 



FIGURE 2.1 One-dimensional user position. 
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In GPS the position of the satellite is known from the ephemeris data trans- 
mitted by the satellite. One can measure the distance from the receiver to the 
satellite. Therefore, the position of the receiver can be determined. 

In the above discussion, the distance measured from the user to the satellite 
is assumed to be very accurate and there is no bias error. However, the distance 
measured between the receiver and the satellite has a constant unknown bias, 
because the user clock usually is different from the GPS clock. In order to 
resolve this bias error one more satellite is required. Therefore, in order to find 
the user position five satellites are needed. 

If one uses four satellites and the measured distance with bias error to mea- 
sure a user position, two possible solutions can be obtained. Theoretically, one 
cannot determine the user position. However, one of the solutions is close to the 
earth’s surface and the other one is in space. Since the user position is usually 
close to the surface of the earth, it can be uniquely determined. Therefore, the 
general statement is that four satellites can be used to determine a user position, 
even though the distance measured has a bias error. 

The method of solving the user position discussed in Sections 2.5 and 2.6 
is through iteration. The initial position is often selected at the center of the 
earth. The iteration method will converge on the correct solution rather than 
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the one in space. In the following discussion four satellites are considered the 
minimum number required in finding the user position. 



2.4 BASIC EQUATIONS FOR FINDING USER POSITION 

In this section the basic equations for determining the user position will be pre- 
sented. Assume that the distance measured is accurate and under this condition 
three satellites are sufficient. In Figure 2 . 3 , there are three known points at loca- 
tions ri or (x\ , vi, zi), r 3 or (%2, yi, Zi), and r 3 or (*3, y 3 , Z3), and an unknown 
point at r u or (x u , y u , z u )- If the distances between the three known points to 
the unknown point can be measured as p\, p2, and pj, these distances can be 
written as 



Pi = v 7 (x! - x u ) 2 + (yi - y u ) 2 + (zi - Zu) 2 

P2 = V (x 2 - x u ) 2 + (y 2 ~ y u ) 2 + (Z2 - Zu) 2 

P3 = V (x 3 - x u ) 2 + (y 3 - y u ) 2 + (Z3 ~ Zu) 2 ( 2 . 1 ) 

Because there are three unknowns and three equations, the values of x u , y u , 

and z u can be determined from these equations. Theoretically, there should be 




FIGURE 2.3 Use three known positions to find one unknown position. 
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two sets of solutions as they are second-order equations. Since these equations 
are nonlinear, they are difficult to solve directly. However, they can be solved 
relatively easily with linearization and an iterative approach. The solution of 
these equations will be discussed later in Section 2.6. 

In GPS operation, the positions of the satellites are given. This information 
can be obtained from the data transmitted from the satellites and will be dis- 
cussed in Chapter 5. The distances from the user (the unknown position) to 
the satellites must be measured simultaneously at a certain time instance. Each 
satellite transmits a signal with a time reference associated with it. By measur- 
ing the time of the signal traveling from the satellite to the user the distance 
between the user and the satellite can be found. The distance measurement is 
discussed in the next section. 



2.5 MEASUREMENT OF PSEUDORANGE< 2 ) 

Every satellite sends a signal at a certain time t sl . The receiver will receive the 
signal at a later time t u . The distance between the user and the satellite i is 

PiT = c(t u -t si ) (2.2) 

where c is the speed of light, prr is often referred to as the true value of pseu- 
dorange from user to satellite i, t S i is referred to as the true time of transmission 
from satellite i, t u is the true time of reception. 

From a practical point of view it is difficult, if not impossible, to obtain the 
correct time from the satellite or the user. The actual satellite clock time t' si and 
actual user clock time t' u are related to the true time as 

t' si = t si + Abi 

C -t u + b ut (2.3) 

where A b t is the satellite clock error, b ut is the user clock bias error. Besides 
the clock error, there are other factors affecting the pseudorange measurement. 
The measured pseudorange p, can be written as (2) 

Pi= PiT+ A D, - c(Abi - b ut ) + c(ATt + A/, + u, + Au ,) (2.4) 

where AD, is the satellite position error effect on range, AT, is the tropospheric 
delay error, A 7, is the ionospheric delay error, a, is the receiver measurement 
noise error, Av t is the relativistic time correction. 

Some of these errors can be corrected; for example, the tropospheric delay 
can be modeled and the ionospheric error can be corrected in a two-frequency 
receiver. The errors will cause inaccuracy of the user position. However, the 
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user clock error cannot be corrected through received information. Thus, it will 
remain as an unknown. As a result, Equation (2.1) must be modified as 



Pi =V(xi- x u ) 2 + ( yi - y u ) 2 + (zi - Zu) 2 + b u 

Pi = V (x 2 - x u ) 2 + ( y 2 - y u ) 2 + (22 - Zu) 2 + b u 

P 3 = v 7 (x 3 - x u ) 2 + ( y 3 - y u ) 2 + (Z 3 - Zu) 2 + b u 

P 4 = V 7 (x 4 - x u ) 2 + ( y 4 - y u ) 2 + (Z4 ~ Zu) 2 + b u (2.5) 

where b u is the user clock bias error expressed in distance, which is related to 
the quantity b ut by b u — cb ut . In Equation (2.5), four equations are needed to 
solve for four unknowns x u , y u , z u , and b u . Thus, in a GPS receiver, a min- 
imum of four satellites is required to solve for the user position. The actual 
measurement of the pseudorange will be discussed in Chapter 9. 



2.6 SOLUTION OF USER POSITION FROM PSEUDORANGES 

It is difficult to solve for the four unknowns in Equation (2.5), because they 
are nonlinear simultaneous equations. One common way to solve the problem 
is to linearize them. The above equations can be written in a simplified form 
as 



Pi = V( Xi - x u ) 2 + ( yi - y u ) 2 + (zi - Zu) 2 + b u (2.6) 

where i = 1, 2, 3, and 4, and x u , y u , z u , and b u are the unknowns. The pseudo- 
range pi and the positions of the satellites x,, y ; , zi are known. 

Differentiate this equation, and the result is 

(xi - x u )bx u + ( y t - y u )hu + ( Zi - z u )bz u , , 

opt , E ob u 

Vixi-Xu^ + iyi-yuf + Ui-Zu) 2 



(xi-x u )8x u + (yi-y u )8y u + (zi-Zu)8z u , s , „ 

= 7 + ob u (2.7) 

Pi ~ b u 

In this equation, 8x u , 8y u , 8z u , and 8b u can be considered as the only unknowns. 
The quantities x u , y u , z u , and b u are treated as known values because one can 
assume some initial values for these quantities. From these initial values a new 
set of 8x u , 8y u , 8z u , and 8b u can be calculated. These values are used to modify 
the original x u , y u , z u , and b u to find another new set of solutions. This new set 
of x u , y u , z u , and b u can be considered again as known quantities. This process 
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continues until the absolute values of bx u , 8y u , bz u , and bb u are very small and 
within a certain predetermined limit. The final values of x u ,y u ,z u , and b u are 
the desired solution. This method is often referred to as the iteration method. 

With bx u , by u , bz u , and bb u as unknowns, the above equation becomes a 
set of linear equations. This procedure is often referred to as linearization. The 
above equation can be written in matrix form as 



~bpi - 




" «n ai2 ai3 1 


~ bx u 


bp 2 




<*21 <*22 <*23 1 


by u 


<$P3 




CC31 CC32 CX33 1 


&Zu 


-bpA- 




_ CX41 «42 0!43 1 - 


-bb u 



(2.8) 



where 



a i2 



y< - y» 

Pi ~ K 



<*i3 



Zj - Zu 
Pi - K 



(2.9) 



The solution of Equation (2.8) is 



~ bx u ~ 




"an ai2 ai 3 1 


- 1 




by u 




<*21 <*22 <*23 1 




bp 2 


bz u 




a3i a32 a33 1 




<5/03 


-bb u _ 




_ a4i a42 a43 1 _ 




-<5,04- 



where [ ] 1 represents the inverse of the a matrix. This equation obviously does 
not provide the needed solutions directly; however, the desired solutions can be 
obtained from it. In order to find the desired position solution, this equation must 
be used repetitively in an iterative way. A quantity is often used to determine 
whether the desired result is reached and this quantity can be defined as 



bv - 



bxl + byl + bzl + bb 2 u 



( 2 . 11 ) 



When this value is less than a certain predetermined threshold, the iteration will 
stop. Sometimes, the clock bias b u is not included in Equation (2.11). 

The detailed steps to solve the user position will be presented in the next 
section. In general, a GPS receiver can receive signals from more than four 
satellites. The solution will include such cases as when signals from more than 
four satellites are obtained. 
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2.7 POSITION SOLUTION WITH MORE THAN FOUR SATELLITES^) 

When more than four satellites are available, a more popular approach to solve 
the user position is to use all the satellites. The position solution can be obtained 
in a similar way. If there are n satellites available where n > 4, Equation (2.6) 
can be written as 



Pi = v 7 (xt - x u ) 2 + ( y, - y u ) 2 + (zi - z u ) 2 + b u (2.12) 



where i = 1, 2, 3, . . . n. The only difference between this equation and Equation 
(2.6) is that n > 4. 

Linearize this equation, and the result is 



~ dpi ~ 

dp 2 

&P3 

dp4 




"an an ai 3 1" 

«21 0122 0123 1 

a 3 i a 32 a 33 1 

a 4 i a 42 a 43 1 




" dx u ~ 

by u 

bz u 


— bpn — 




_a„i a n2 a „3 1_ 




_ db u _ 



where 



(2.13) 



a n = 



Xi - x u 
Pi ~ b u 



yi - y u 

Pi ~ b u 



Pi ~ b u 



(2.9) 



Equation (2.13) can be written in a simplified form as 

dp-a8x (2.14) 

where 8p and 8x are vectors, a is a matrix. They can be written as 

dp = [5pi dp2 ■ ■ ■ dp n ] T 
8x = [8x u 8y u bz u &b u ] T 

ail «12 «13 1 

«21 0122 «23 1 

«31 «32 a33 1 . 

a 4 l «42 «43 1 I ( 2 - 15 ) 

L a«i a«2 a«3 1 . 
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where [ ] r represents the transpose of a matrix. Since a is not a square matrix, 
it cannot be inverted directly. Equation (2.13) is still a linear equation. If there 
are more equations than unknowns in a set of linear equations, the least-squares 
approach can be used to find the solutions. The pseudoinverse of the a can be 
used to obtain the solution. The solution is® 



8x - [a T a\ l a r 8p 



(2.16) 



From this equation, the values of 8x u , 8y u , 8z u , and 8b u can be found. In general, 
the least-squares approach produces a better solution than the position obtained 
from only four satellites, because more data are used. 

The following steps summarize the above approach: 

A. Choose a nominal position and user clock bias x u q, y u o, z u o> b u o to rep- 
resent the initial condition. For example, the position can be the center 
of the earth and the clock bias zero. In other words, all initial values are 
set to zero. 

B. Use Equation (2.5) or (2.6) to calculate the pseudorange p,. These p, val- 
ues will be different from the measured values. The difference between 
the measured values and the calculated values is <5p,. 

C. Use the calculated p, in Equation (2.9) to calculate an, an, an. 

D. Use Equation (2.16) to find 8x u , 8y u , 8z u , 8b u . 

E. From the absolute values of 8x u , 8y u , 8z u , 8b u and from Equation (2.11) 
calculate 8v. 

F. Compare 8v with an arbitrarily chosen threshold; if 8v is greater than the 
threshold, the following steps will be needed. 

G. Add these values 8x u , 8y u , 8z u , 8b u to the initial chosen position x u o, 
y u o, ZuO, and the clock bias b u o; a new set of positions and clock bias 
can be obtained and they will be expressed as x u \, y u \ , z u i, b u \. These 
values will be used as the initial position and clock bias in the following 
calculations. 

H. Repeat the procedure from A to G, until 8v is less than the threshold. The 
final solution can be considered as the desired user position and clock 
bias, which can be expressed as x u , y u , z u , b u . 



In general, the 8v calculated in the above iteration method will keep decreas- 
ing rapidly. Depending on the chosen threshold, the iteration method usually can 
achieve the desired goal in less than 10 iterations. A computer program (p2_l) 
to calculate the user position is listed at the end of this chapter. 
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2.8 USER POSITION IN SPHERICAL COORDINATE SYSTEM 

The user position calculated from the above discussion is in a Cartesian coor- 
dinate system. It is usually desirable to convert to a spherical system and label 
the position in latitude, longitude, and altitude as the every-day maps use these 
notations. The latitude of the earth is from -90 to 90 degrees with the equator 
at 0 degree. The longitude is from -180 to 180 degrees with the Greenwich 
meridian at 0 degree. The altitude is the height above the earth’s surface. If 
the earth is a perfect sphere, the user position can be found easily as shown 
in Figure 2.4. From this figure, the distance from the center of the earth to the 
user is 



r - ‘\J x u + yi + 



(2.17) 



The latitude L c is 



L c - tan 1 



(2.18) 



The longitude l is 
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l - tan 



(2.19) 



The altitude h is 



( 2 . 20 ) 



where r e is the radius of an ideal spherical earth or the average radius of the 
earth. Since the earth is not a perfect sphere, some of these equations need to 
be modified. 



2.9 EARTH GEOMETRY! 4-6 ) 

The earth is not a perfect sphere but is an ellipsoid; thus, the latitude and altitude 
calculated from Equations (2.18) and (2.20) must be modified. However, the 
longitude / calculated from Equation (2.19) also applies to the nonspherical 
earth. Therefore, this quantity does not need modification. Approximations will 
be used in the following discussion, which is based on references 4 through 6. 
For an ellipsoid, there are two latitudes. One is referred to as the geocentric 
latitude L c , which is calculated from the previous section. The other one is the 
geodetic latitude L and is the one often used in every-day maps. Therefore, the 
geocentric latitude must be converted to the geodetic latitude. Figure 2.5 shows 
a cross section of the earth. In this figure the x-axis is along the equator, the 
y-axis is pointing inward to the paper, and the z-axis is along the north pole of 
the earth. Assume that the user position is on the x-z plane and this assumption 
does not lose generality. The geocentric latitude L c is obtained by drawing a 
line from the user to the center of the earth, which is calculated from Equation 
(2.18). 

The geodetic latitude is obtained by drawing a line perpendicular to the sur- 
face of the earth that does not pass the center of the earth. The angle between 
this line and the x is the geodetic latitude L. The height of the user is the dis- 
tance h perpendicular and above the surface of the earth. 

The following discussion is used to determine three unknown quantities from 
two known quantities. As shown in Figure 2.5, the two known quantities are 
the distance r and the geocentric latitude L c and they are measured from the 
ideal spherical earth. The three unknown quantities are the geodetic latitude 
L, the distance r 0 , and the height h. All three quantities are calculated from 
approximation methods. Before the actual calculations of the unknowns, let us 
introduce some basic relationships in an ellipse. 
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2.10 BASIC RELATIONSHIPS IN AN ELLIPSE* 4 ' 7 ) 

In order to derive the relationships mentioned in the previous section, it is con- 
venient to review the basic functions in an ellipse. Figure 2.6 shows an ellipse 
which can be used to represent a cross section of the earth passing through the 
polar axis. 

Let us assume that the semi-major axis is a e , the semi-minor axis is b e , and 
the foci are separated by 2 c e . The equation of the ellipse is 



a\ b] 



= 1 and 



( 2 . 21 ) 



The eccentricity e e is defined as 



The ellipticity e p is defined as 






(2.23) 



where a e = 6378137 ± 2 m , b e = 6356752.3142 m, e e = 0.0818191908426, and 
e p - 0.0033528 1 066474. (6 - 7 * The value of b e is calculated from a e ; thus, the 
result has more decimal points. 

From the user position P draw a line perpendicular to the ellipse that inter- 
cepts it at A and the x-axis at C. To help illustrate the following relation a circle 
with radius equal to the semi-major axis a e is drawn as shown in Figure 2.6. 
A line is drawn from point A perpendicular to the x-axis and intercepts it at E 
and the circle at D. The position A(x, y) can be found as 

x = OE - OD cos (8 = a e cos /3 

Z=AE = DE — ={a e sin p) — =b e sin j8 (2.24) 

Q-e 

The second equation can be obtained easily from the equation of a circle x 2 + 
y 2 - a 2 and Equation (2.21). The tangent to the ellipse at A is dz/dx. Since line 
CP is perpendicular to the tangent, 
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(2.25) 



From these relations let us find the relation between angle 0 and L. Taking the 
derivative of x and z of Equation (2.24), the results are 

dx - -a e sin j3d(3 

dz - b e cos I3d(3 (2.26) 



Thus 



a e 

~b~e 



tan (3 




(2.27) 



From these relationships let us find the three unknowns. 



2.11 CALCULATION OF ALTITUDE^) 

In the following three sections the discussion is based on reference 5. From 
Figure 2.7 the height h can be found from the law of cosine through the triangle 
OPA as 




FIGURE 2.7 Altitude and latitude illustration. 
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r 2 =r q - 2 r 0 h cos(7 r - D 0 ) + /i 2 = r 2 + 2r 0 /i cos D 0 + /i 2 (2.28) 

where is the distance from the center of the earth to the point on the surface 
of the earth under the user position. The amplitude of r can be found from 
completing the square for r$ + h and taking the square root as 



r = [(r 0 +/i) 2 -2r 0 /i(l-cosD 0 )] 1/2 = (r 0 +ft)|l - + ^ ^ ] ( 2 - 29 ) 

Since angle Do is very small, it can be approximated as 
Dl 

1 - cos D 0 = (2.30) 
where Dq is the angle expressed in radians. The r value can be written as 



r « (to + h) 



2hr 0 D 2 /2 

( r 0 + h ) 2 



1/2 



ro + h 



hr 0 D 2 
2 (r 0 + h) 



(2.31) 



At latitude of 45 degrees Do (« 1/297 radian) becomes maximum. If Do is 
neglected, the result is 



r « r 0 + ft - r ° 0 “ r 0 + h (2.32) 

2(r 0 + h) 

Using this result, if h = 100 km, and t'Q - r f , - 6368 km (the average radius of 
the earth), the error term calculated is less than 0.6 m. Thus 

h = r-r 0 (2.33) 

is a good approximation. However, in this equation the value of r 0 must be 
evaluated, as discussed in Section 2.12. 



2.12 CALCULATION OF GEODETIC LATITUDE* 5 - 7 ) 

Referring to Figure 2.7, the relation between angles L and L c can be found from 
the triangle OPC. From the simple geometry it can be seen that 



L=L C + D 



(2.34) 
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If the angle D can be found, the relation between L and L c can be obtained. To 
find this angle, let us find the distance OC first. Combining Equations (2.24) 
and (2.27), the following result is obtained: 



AE b e sin |8 

= a e cos p 

tan L tan L 



OC = OE- CE = a e 

- a e cos /3[1 - (1 - e 2 )\ - a e e 2 e cos j3 = e 2 e OE 



where /3 is not shown in this figure but is shown in Figure 2.6. 
From the triangle OPC and the law of sine, one can write 



From Equation (2.35), 



but 



sinD sin(7r-L) 
OC = r 



OC = e~OE = e 2 ro cos L co 



L co =L- D 0 

Therefore, 

OC = e 2 ro cos (L - D a ) = e^ro(cos L cos Dq + sin L sin Do) 
From Equation (2.23), the ellipticity e p of the earth is 
a e - b e 



(2.39) 



The eccentricity and the ellipticity can be related as 



a 2 - b 2 e (a e - b e ) (a e + b e ) 



(2 a e -a e + b e ) 



= e p (2-e p ) (2.41) 



Substituting Equations (2.39) and (2.41) into Equation (2.36), the result is 



sin D = 2e„ ( 1 — ] — — sin 2 L cos D 0 + sin 2 L sin D 0 ) 

y \ 2 J r 0 + h \ 2 J 



(2.42) 
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In the above equation the relation r = r 0 + h is used. Since D and D 0 are both 
very small angles, the above equation can be written as 

D = 2cp(l--y) sin 2L + D 0 sin 2 l)j (2.43) 

The relations 

sinD ~ D; sin Do = Do cos Do ~ 1 (2.44) 

are used in obtaining the results of Equation (2.43). If the height h = 0, then 
from Figure 2.7 D = Dq. Using this relation Equation (2.43) can be written as 



Dq = e„ sin 2L + e 



(2.45) 



where 



e i = - sin L + 2e 2 sin 2 L sin 2 L + ... < 1 .6 arc - sec (2.46) 

Substitute the approximation of Dq ~ e p sin 2 L into Equation (2.43) and the 
result is 



D = 2e p ^ 1 - -y) ( 1 - y-j ^ j sin 2L + e p sin 2L sin 2 L j (2.47) 

or 

D = e p sin 2L + e (2.48) 

where 



e = - — sin 2L - — sin 2L H (2.49) 

2 r 0 

This error is less than 4.5 arc-sec for h = 30 km. Using the approximate value 
of D, the relation between angle L and L c can be found from Equation (2.34) 
as 
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L = L C + e p sin 2L (2.50) 

This is a nonlinear equation that can be solved through the iteration method. 
This equation can be written in a slightly different form as 

L i+i - L c + e p sin2Lj (2.51) 

where i - 0, 1,2,... . One can start with Lq - L c . If the difference (L,- + 1 - L,) is 
smaller than a predetermined threshold, the last value of L, can be considered 
as the desired one. It should be noted that during the iteration method L c is a 
constant that is obtained from Equation (2.18). 



2.13 CALCULATION OF A POINT ON THE SURFACE OF THE EARTH( 5 ) 

The final step of this calculation is to find the value ro in Equation (2.33). This 
value is also shown in Figure 2.7. The point A ( x , y) is on the ellipse; therefore, 
it satisfies the following elliptic Equation (2.21). This equation is rewritten here 
for convenience, 



b 2 



(2.52) 



where a e and b e are the semi-major and semi-minor axes of the earth. From 
Figure 2.7, the x and y values can be written as 

x = ro cos L co 

y = r 0 sin L co (2.53) 

Substituting these relations into Equation (2.52) and solving for ro, the result 
is 



sin 2 L co ' 


\ 2 ( bl COS 2 L co ■ 


b] , 


1 °[ 




a 2 e b 2 e 


a 2 e 1- 


( 1 -^) cosi “ 



r 0 = b e M + — e 2 e cos 2 L co + ■ ■ ■ 




SELECTION 



Use Equation (2.23) to replace b e by a e , Equation (2.41) to replace e e by e p , 
and L to replace L co because L~ L co , and then 




« a e ( 1 - e p ){\ +e p -e p sin 2 L + • • •) (2.55) 

In this equation the higher order of e p is neglected. The value of rg can be 
found as 

r 0 ~ a e ( 1 - e p sin 2 L ) (2.56) 

To solve for the latitude and altitude of the user, use Equation (2.51) to find 
the geodetic latitude L first. Then use Equation (2.56) to find rg, and finally, 
use Equation (2.33) to find the altitude. The result is 

h~r - r 0 ~ yjxl+yl +zl - a e (l - e p sin 2 L) (2.57) 



2.14 SATELLITE SELECTION^ 1 ' 8 ) 

A GPS receiver can simultaneously receive signals from 4 up to 11 satellites, 
if the receiver is on the surface of the earth. Under this condition, there are 
two approaches to solve the problem. The first one is to use all the satellites to 
calculate the user position. The other approach is to choose only four satellites 
from the constellation. The usual way is to utilize all the satellites to calculate 
the user position, because additional measurements are used. In this section and 
section 2.15 the selection of satellites will be presented. In order to focus on 
this subject only the four-satellite case will be considered. 

If there are more than four satellite signals that can be received by a GPS 
receiver, a simple way is to choose only four satellites and utilize them to solve 
for the user position. Under this condition, the question is how to select the four 
satellites. Let us use a two-dimensional case to illustrate the situation, because 
it is easier to show graphically. In order to solve a position in a two-dimen- 
sional case, three satellites are required considering the user clock bias. In this 
discussion, it is assumed that the user position can be uniquely determined as 
discussed in Section 2.3. If this assumption cannot be used, four satellites are 
required. 

Figure 2.8a shows the results measured by three satellites on a straight line, 
and the user is also on this line. Figure 2.8b shows that the three satellites 
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and the user form a quadrangle. Two circles with the same center but different 
radii are drawn. The solid circle represents the distance measured from the user 
to the satellite with bias clock error. The dotted circle represents the distance 
after the clock error correction. From observation, the position error in Figure 
2.8a is greater than that in Figure 2.8b because in Figure 2.8a all three dotted 
circles are tangential to each other. It is difficult to measure the tangential point 
accurately. In Figure 2.8b, the three circles intersect each other and the point 
of intersection can be measured more accurately. Another way to view this 
problem is to measure the area of a triangle made by the three satellites. In 
Figure 2.8a the total area is close to zero, while in Figure 2.8b the total area is 
quite large. In general, the larger the triangle area made by the three satellites, 
the better the user position can be solved. 

The general rule can be extended to select the four satellites in a three-dimen- 
sional case. It is desirable to maximize the volume defined by the four satellites. 
A tetrahedron with an equilateral base contains the maximum volume and there- 
fore can be considered as the best selection. Under this condition, one satellite 
is at zenith and the other three are close to the horizon and separated by 120 
degrees.' ® This geometry will generate the best user position estimation. If all 
four satellites are close to the horizon, the volume defined by these satellites 
and the user is very small. Occasionally, the user position error calculated with 
this arrangement can be extremely large. In other words, the 8v calculated from 
Equation (2.11) may not converge. 



2.15 DILUTION OF PRECISION^) 

The dilution of precision (DOP) is often used to measure user position accuracy. 
There are several different definitions of the DOP. All the different DOPs are 
a function of satellite geometry only. The positions of the satellites determine 
the DOP value. A detailed discussion can be found in reference 8. Here only 
the definitions and the limits of the values will be presented. 

The geometrical dilution of precision (GDOP) is defined as 

GDOP = yja} + <Jy + a 2 + a\ (2.58) 

where a is the measured rms error of the pseudorange, which has a zero mean, 
o x o y o z are the measured rms errors of the user position in the xyz directions, 
and Ob is the measured rms user clock error expressed in distance. 

The position dilution of precision is defined as 



PDOP - ^ a l + Oy + a z 2 



(2.59) 
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The horizontal dilution of precision is defined as 

HDOP = — A /ff„ 2 + a v 2 
a V ' 



(2.60) 



The vertical dilution of precision is 



VDOP = — (2.61) 

<7 

The time dilution of precision is 

TDOP = — (2.62) 

a 

The smallest DOP value means the best satellite geometry for calculating 
user position. It is proved in reference 8 that in order to minimize the GDOP, 
the volume contained by the four satellites must be maximized. Assume that 
the four satellites form the optimum constellation. Under this condition the ele- 
vation angle is 0 degree and three of the four satellites form an equilateral tri- 
angle. The observer is at the center of the base of the tetrahedron. Under this 
condition, the DOP values are: GDOP = v/3 = 1.73, PDOP = 2\/2/3 ~ 1.63, 
HDOP = VDOP = 2/V3 = 1.15, and TDOP = 1/V3 « 0.58. These values can 
be considered as the minimum values (or the limits) of the DOPs. In selecting 
satellites, the DOP values should be as small as possible in order to generate 
the best user position accuracy. 



2.16 SUMMARY 

This chapter discusses the basic concept of solving the GPS user position. First 
use four or more satellites to solve the user position in terms of latitude, lon- 
gitude, altitude, and the user clock bias as discussed in Section 2.5. However, 
the solutions obtained through this approach are for a spherical earth. Since 
the earth is not a perfect sphere, the latitude and altitude must be modified to 
reflect the ellipsoidal shape of the earth. Equations (2.51) and (2.57) are used 
to derive the desired values. These results are shown in Figure 2.9 as a quick 
reference. Finally, the selection of satellites and the DOP are discussed. 
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FIGURE 2.9 Relations to change from spherical to ellipsoidal earth. 
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% p2_l .m 

% Userpos.m use pseudorange and satellite positions to calculate user 

position 

% JT 30 April 96 

% ***** input data ***** 

sp (1 : 3 , l:nsat) ; % satellite position which has the following format 



r: 



Y% 

Yi 

Yt, 



% 

Z2 

Z3 



pr (l:nsat) ; % is the measured pseudo-range which has the format as 
% pr=[prl pr2 pr3 . . . prnn] T ; 

nn=nsat; % is the number of satellites 

% ***** select initial guessed positions and clock bias ***** 
x_guess = 0; y_guess = 0; z_guess = 0; bu = 0; 
gu(l) =x_guess; gu(2) =y_guess; gu(3) = z guess; 

% Calculating rao the pseudo-range as shown in Equation (2.1) the 
% clock bias is not included 
for j = 1 :nsat 

rao(j)=( (gu (1 ) -sp(l, j) ) A 2+(gu(2) -sp(2, j}) A 2+ (gu(3) 

-sp (3 , j ) ) A 2 ) A . 5 ; 
end 

% generate the fourth column of the alpha matrix in Eq. 2.15 
alpha!:, 4) = ones (nsat, 1) ; 

erro=l ; 

while erroi . 01; 
for j = 1 :nsat; 
f or k = 1 : 3 ; 

alpha (j,k) = ( gu ( k ) - sp ( k , j ) ) / ( rao ( ^ ) J % find first 
%3 columns of alpha matrix 
end 
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end 

drao = pr - (rao + ones (l,nsat) *bu) ;%** find delta rao 

% includes clock bias 
dl = pinv(h) *drao' ; % Equation (2.16) 
bu = bu + dl (4) ; % new clock bias 
f or k = 1 : 3 ; 

gu(k) =gu(k) +dl(k); %**find new position 
end 

erro=dl (1) A 2+dl (2) A 2+dl (3) A 2; % find error 
for j = 1 :nsat; 

rao ( j ) = ( ( gu ( 1 ) -sp (1, j ) ) A 2+ (gu (2) -sp (2 , j ) ) A 2+ (gu (3 ) - 
sp (3 , j ) ) A 2) A . 5 ; % find new rao without clock bias 

end 

end 

% ***** Final result in spherical coordinate system ***** 

xuser = gu(l) ; yuser = gu(2) ; zuser = gu(3) ; bias = bu; 

rsp = (xuser A 2+yuser A 2+zuser A 2 ) a . 5; % radius of spherical earth 
% Eq 2.17 

Lc = atan( zuser/ (xuser A 2+yuser A 2) A .5) ; % latitude of spherical 
% earth Eq 2 . 18 

lsp = atan(yuser/xuser) *180/pi; % longitude spherical and flat 
% earth Eq 2.19 

% ***** Converting to practical earth shape ***** 

e=l/298. 257223563; 

Ltemp=Lc ; 
errol=l ; 

while errol>le-6; % calculating latitude by Eq. 2.51 
L=Lc+e*sin(2*Ltemp) ; 
errol=abs (Ltemp-L) ; 

Ltemp^L; 

end 

Lflp=L* 180/pi ; % latitude of flat earth 
re=6378137 ; 

h=rsp-re* (1-e* (sin (L) A 2 ) ) ; % altitude of flat earth 

lsp = lsp; % longitude of flat earth 

upos = [xuser yuser zuser bias rsp Lflp lsp h] ' ; 
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CHAPTER. THREE 



Satellite Constellation 



3.1 INTRODUCTION 

The previous chapter assumes that the positions of the satellites are known. 
This chapter will discuss the satellite constellation and the determination of the 
satellite positions. Some special terms related to the orbital mechanics, such as 
sidereal day, will be introduced. The satellite motion will have an impact on 
the processing of the signals at the receiver. For example, the input frequency 
shifts as a result of the Doppler effect. Such details are important for the design 
of acquisition and tracking loops in the receiver. However, in order to obtain 
some of this information a very accurate calculation of the satellite motion is not 
needed. For example, the actual orbit of the satellite is elliptical but it is close 
to a circle. The information obtained from a circular orbit will be good enough 
to find an estimation of the Doppler frequency. Based on this assumption the 
circular orbit is used to calculate the Doppler frequency, the rate of change of 
the Doppler frequency, and the differential power level. From the geometry of 
the satellite distribution, the power level at the receiver can also be estimated 
from the transmission power. This subject is presented in the final section in this 
chapter. 

In order to find the location of the satellite accurately, a circular orbit is insuf- 
ficient. The actual elliptical satellite orbit must be used. Therefore, the complete 
elliptical satellite orbit and Kepler’s law will be discussed. Information obtained 
from the satellite through the GPS receiver via broadcast navigation data such 
as the mean anomaly does not provide the location of the satellite directly. 
However, this information can be used to calculate the precise location of the 
satellite. The calculation of the satellite position from these data will be dis- 
cussed in detail. 
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3.2 CONTROL SEGMENT OF THE GPS SYSTEM^ -3 ) 

This section will provide a very brief idea of the GPS system. The GPS sys- 
tem may be considered as comprising three segments: the control segment, the 
space segment, and the user segment. The space segment contains all the satel- 
lites, which will be discussed in Chapters 3, 4, and 5. The user segment can 
be considered the base of receivers and their processing, which is the focus of 
this text. The control segment will be discussed in this section. 

The control segment consists of five control stations, including a master con- 
trol station. These control stations are widely separated in longitude around the 
earth. The master control station is located at Falcon Air Force Base, Colorado 
Springs, CO. Operations are maintained at all times year round. The main pur- 
pose of the control stations is to monitor the performance of the GPS satellites. 
The data collected from the satellites by the control stations will be sent to 
the master control station for processing. The master control station is respon- 
sible for all aspects of constellation control and command. A few of the oper- 
ation objectives are presented here: (1) Monitor GPS performance in support 
of all performance standards. (2) Generate and upload the navigation data to 
the satellites to sustain performance standards. (3) Promptly detect and respond 
to satellite failure to minimize the impact. Detailed information on the control 
segment can be found in reference 3. 

3.3 SATELLITE CONSTELLATION^ -9 ) 

There are a total of 24 GPS satellites divided into six orbits and each orbit has four 
satellites. Each orbit makes a 55-degree angle with the equator, which is referred 
to as the inclination angle. The orbits are separated by 60 degrees to cover the 
complete 360 degrees. The radius of the satellite orbit is 26,560 km and it rotates 
around the earth twice in a sidereal day. Table 3.1 lists all these parameters. 

The central body of the Block HR satellite is a cube of approximately 6 ft 
on each side.® The span of the solar panel is about 30 ft. The lift-off weight 
of the spacecraft is 4,480 pounds and the on-orbit weight is 2,370 pounds. 



TABLE 3.1 Characteristics of GPS Satellites 

Constellation 



Number of satellites 
Number of orbital planes 
Number of satellites per orbit 
Orbital inclination 
Orbital radius® 

Period® 

Ground track repeat 



24 

6 

4 

55° 

26560 km 

11 hrs 57 min 57.26 sec 
sidereal day 
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The four satellites in an orbit are not equally spaced. Two satellites are sepa- 
rated by 30.0-32.1 degrees. The other two make three angles with the first two 
satellites and these angles range 92.38-130.98 degrees. (9) The spacing has been 
optimized to minimize the effects of a single satellite failure on system degra- 
dation. At any time and any location on the earth, neglecting obstacles such as 
mountains and tall buildings, a GPS receiver should have a direct line of sight 
and be receiving signals from 4 to 11 satellites. A majority of the time a GPS 
receiver can receive signals from more than four satellites. Since four satellites 
are the minimum required number to find the user position, this arrangement can 
provide user position at any time and any location. For this 24-satellite constel- 
lation with a 5-degree elevation mask angle, more than 80% of the time seven 
or more satellites are in view. (9) A user at 35 degrees latitude corresponds to 
the approximate worst latitude where momentarily there are only four satellites 
in view (approximately .4% of the time). 

The radius of the earth is 6,378 km around the equator and 6,357 km passing 
through the poles, and the average radius can be considered as 6,368 km. The 
radius of the satellite orbit is 26,560 km, which is about 20,192 km (26,560 - 
6,368) above the earth’s surface. This height agrees well with references 6 and 
9. This height is approximately the shortest distance between a user on the sur- 
face of the earth and the satellite, which occurs at around zenith or an elevation 
angle of approximately 90 degrees. Most GPS receivers are designed to receive 
signals from satellites above 5 degrees. For simplicity, let us assume that the 
receiver can receive signals from satellites at the zero-degree point. The distance 
from a satellite on the horizon to the user is 25,785 km (V 26, 560 2 - 6, 368 2 ). 
These distances are shown in Figure 3.1. 

From the distances in Figure 3.1 one can see that the time delays from the 
satellites are in the range of 67 ms (20,192 km/c) to 86 ms (25,785 km/c), 
where c is the speed of light. If the user is on the surface of the earth, the 
maximum differential delay time from two different satellites should be within 
19 (86-67) ms. In this figure, the angle a is approximately 76.13 degrees and 
the angle /3 is approximately 13.87 degrees. Therefore, the transmitting antenna 
on the satellite need only have a solid angle of 13.87 degrees to cover the earth. 
However, the antenna for the L| band is 21.3 degrees and the antenna for the 
L 2 band is 23.4 degrees. Both are wider than the minimum required angle. The 
solid angle of 21.3 degrees will be used in Section 3.13 to estimate the power 
to the receiver. The antenna pattern will be further discussed in Section 5.2. 



3.4 MAXIMUM DIFFERENTIAL POWER LEVEL FROM DIFFERENT 
SATELLITES 

From Figure 3.1 one can calculate the relative power level of the received sig- 
nals on the surface of the earth. The transmitting antenna pattern is designed to 
directly aim at the earth underneath it. However, the distances from the receiver 
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to various satellites are different. The shortest distance to a satellite is at zenith 
and the farthest distance to a satellite is at horizon. Suppose the receiver has an 
omnidirectional antenna pattern. Since the power level is inversely proportional 
to the distance square, the difference in power level can be found as 

A P = 101 ° g (|^^2") ~2.\dB (3.1) 

It is desirable to receive signals from different satellites with similar strength. 
In order to achieve this goal, the transmitting antenna pattern must be properly 
designed. The beam is slightly weaker at the center to compensate for the power 
difference. 



3.5 SIDEREAL DAY* 10 ' 11 ) 

Table 3.1 indicates that the satellite rotates around the earth twice in a side- 
real day. The sidereal day is slightly different from an apparent solar day. The 
apparent day has 24 hours and it is the time used daily. The apparent solar 
day is measured by the time between two successive transits of the sun across 
our local meridian, because we use the sun as our reference. A sidereal day is 
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FIGURE 3.2 Configuration of apparent solar day and sidereal day. 



the time for the earth to turn one revolution. Figure 3.2 shows the difference 
between the apparent solar day and a sidereal day. In this figure, the effect is 
exaggerated and it is obvious that a sidereal day is slightly shorter than a solar 
day. The difference should be approximately equal to one day in one year which 
corresponds to about 4 min (24 x 60/365) per day. The mean sidereal day is 
23 hrs, 56 min, 4.09 sec. The difference from an apparent day is 3 min, 55.91 
sec. Half a sidereal day is 11 hrs, 58 min, 2.05 sec. This is the time for the 
satellite to rotate once around the earth. From this arrangement one can see 
that from one day to the next a certain satellite will be at approximately the 
same position at the same time. The location of the satellite will be presented 
in the next section. 



3.6 DOPPLER FREQUENCY SHIFT 

In this section, the Doppler frequency shift caused by the satellite motion both 
on the carrier frequency and on the coarse/acquisition (C/A) code will be dis- 
cussed. This information is important for performing both the acquisition and 
the tracking of the GPS signal. 

The angular velocity dd/dt and the speed v s of the satellite can be calculated 
from the approximate radius of the satellite orbit as 
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d9_ 
~di ~ 

V s = 



27 r 

11 x 3600 + 58 x 60 + 2.05 
r s dd 

-V- “26560 km x 1.458 > 
dt 



U 1.458 x 10 4 rad/s 
: 10~ 4 « 3874 m/s 



(3.2) 



where r s is the average radius of the satellite orbit. In 3 min, 55.91 sec, the 
time difference between an apparent solar day and the sidereal day, the satellite 
will travel approximately 914 km (3,874 m/s x 235.91 sec). Referenced to the 
surface of the earth with the satellite in the zenith direction, the corresponding 
angle is approximately .045 radian (914/20.192) or 2.6 degrees. If the satellite 
is close to the horizon, the corresponding angle is approximately .035 radian 
or 2 degrees. Therefore, one can consider that the satellite position changes 
about 2-2.6 degrees per day at the same time with respect to a fixed point 
on the surface of the earth. In Figure 3.3, the satellite is at position S and the 
user is at position A. The Doppler frequency is caused by the satellite velocity 
component Vj toward the user where 



(3.3) 




FIGURE 3.3 Doppler frequency caused by satellite motion. 
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Doppler velocity vs angle 




angle 9 in radian 



FIGURE 3.4 Velocity component toward the user versus angle 0. 



Now let us find this velocity in terms of angle 9. Using the law of cosine in 
triangle OAS, the result is 

AS 2 = r 2 + r 2 - 2 r e r s cos a = r 2 + r 2 - 2 r e r s sin 9 (3.4) 

because of a + 6 = tt/ 2. In the same triangle, using the law of sine, the result 



sin j3 sin /3 r e 

sin a cos 6 AS 

Substituting these results into Equation (3.3), one obtains 

v s r e cos 9 v s r e cos 9 

AS -Vr% + r} - 2 r e r s sin 9 



(3.5) 



(3.6) 



This velocity can be plotted as a function 9 and is shown in Figure 3.4. 

As expected, when 9 = tt/ 2, the Doppler velocity is zero. The maximum 
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Doppler velocity can be found by taking the derivative of v d with respect to 8 
and setting the result equal to zero. The result is 



dv d vr e [r e r s sin 2 8 - (r 2 + r 2 ) sin 8 + r e r s ] 
dd (rj + r 2 - 2 r e r s sin 0) 3 / 2 



Thus sin 8 can be solved 



sin# 



or 8 - sin 



= 0.242 rad 



(3.8) 



At this angle 8 the satellite is at the horizontal position referenced to the 
user. Intuitively, one expects that the maximum Doppler velocity occurs when 
the satellite is at the horizon position and this calculation confirms it. From 
the orbit speed, one can calculate the maximum Doppler velocity iQ m , which 
is along the horizontal direction as 



V s r e 
Vdm = 



3874 x 6368 
26560 



= 929 m/s ~ 2078 miles/h 



(3.9) 



This speed is equivalent to a high-speed military aircraft. The Doppler fre- 
quency shift caused by a land vehicle is often very small, even if the motion 
is directly toward the satellite to produce the highest Doppler effect. For the 
Li frequency (/ = 1575.42 MHz), which is modulated by the C/A signal, the 
maximum Doppler frequency shift is 



fdr- 



frVdm 



1575.42 x 929 
3 x 10 8 



« 4.9 KHz 



(3.10) 



where c is the speed of light. Therefore, for a stationary observer, the maximum 
Doppler frequency shift is around ±5 KHz. 

If a vehicle carrying a GPS receiver moves at a high speed, the Doppler 
effect must be taken into consideration. To create a Doppler frequency shift of 
±5 KHz by the vehicle alone, the vehicle must move toward the satellite at about 
2,078 miles/hr. This speed will include most high-speed aircraft. Therefore, in 
designing a GPS receiver, if the receiver is used for a low-speed vehicle, the 
Doppler shift can be considered as ±5 KHz. If the receiver is used in a high- 
speed vehicle, it is reasonable to assume that the maximum Doppler shift is ±10 
KHz. These values determine the searching frequency range in the acquisition 
program. Both of these values are assumed an ideal oscillator and sampling 
frequency and further discussion is included in Section 6.15. 

The Doppler frequency shift on the C/A code is quite small because of the 
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low frequency of the C/A code. The C/A code has a frequency of 1.023 MHz, 
which is 1,540 (1575.42/1.023) times lower than the carrier frequency. The 
Doppler frequency is 



f c v h 



1.023 x 10 6 x 929 



If the receiver moves at high speed, this value can be doubled to 6.4 Hz. This 
value is important for the tracking method (called block adjustment of syn- 
chronizing signal or BASS program), which will be discussed in Chapter 8. 
In the BASS program, the input data and the locally generated data must be 
closely aligned. The Doppler frequency on the C/A code can cause misalign- 
ment between the input and the locally generated codes. 

If the data is sampled at 5 MHz (referred to as the sampling frequency), each 
sample is separated by 200 ns (referred to as the sampling time). In the tracking 
program it is desirable to align the locally generated signal and the input signal 
within half the sampling time or approximately 100 ns. Larger separation of 
these two signals will lose tracking sensitivity. The chip time of the C/A code 
is 977.5 ns or 1/(1.023 x 10 6 ) sec. It takes 156.3 ms (1/6.4) to shift one cycle 
or 977.5 ns. Therefore, it takes approximately 16 ms (100 x 156.3/977.5) to 
shift 100 ns. In a high-speed aircraft, a selection of a block of the input data 
should be checked about every 16 ms to make sure these data align well with 
the locally generated data. Since there is noise on the signal, using 1 ms of data 
may not determine the alignment accurately. One may extend the adjustment of 
the input data to every 20 ms. For a slow-moving vehicle, the time may extend 
to 40 ms. 

From the above discussion, the adjustment of the input data depends on the 
sampling frequency. Higher sampling frequency will shorten the adjustment 
time because the sampling time is short and it is desirable to align the input 
and the locally generated code within half the sampling time. If the incoming 
signal is strong and tracking sensitivity is not a problem, the adjustment time 
can be extended. However, the input and the locally generated signals should be 
aligned within half a chip time or 488.75 ns (977.5/2). This time can be consid- 
ered as the maximum allowable separation time. With a Doppler frequency of 
6.4 Hz, the adjustment time can be extended to 78.15 ms (1/2 x 6.4). Detailed 
discussion of the tracking program will be presented in Chapter 8. 



3.7 AVERAGE RATE OF CHANGE OF THE DOPPLER FREQUENCY 

In this section the rate of change of the Doppler frequency will be discussed. 
This information is important for the tracking program. If the rate of change 
of the Doppler frequency can be calculated, the frequency update rate in the 
tracking can be predicted. Two approaches are used to find the Doppler fre- 
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quency rate. A very simple way is to estimate the average rate of change of the 
Doppler frequency and the other one is to find the maximum rate of change of 
the Doppler frequency. 

In Figure 3.4, the angle for the Doppler frequency changing from maximum 
to zero is about 1.329 radians (tt/2 -6 - 7t/2 - 0.242). It takes 11 hrs, 58 min, 
2.05 sec for the satellite to travel 2w angle; thus, the time it takes to cover 1.329 
radians is 



t = (11 x 3600 + 58 x 60 + 2.05) =9113 sec (3.12) 

2ir 

During this time the Doppler frequency changes from 4.9 KHz to 0, thus, the 
average rate of change of the Doppler frequency 8f e i r can be simply found as 

4900 

dfdr = nu “ °' 54 Hz/s (3 ' 13) 

This is a very slow rate of change in frequency. From this value a tracking 
program can be updated every few seconds if the frequency accuracy in the 
tracking loop is assumed on the order of 1 Hz. The following section shows 
how to find the maximum frequency rate of change. 



3.8 MAXIMUM RATE OF CHANGE OF THE DOPPLER FREQUENCY 

In the previous section the average rate of change of the Doppler frequency is 
estimated; however, the rate of change is not a constant over that time period. 
In this section we try to find the maximum rate of change of the frequency. The 
rate of change of the speed Vd can be found by taking the derivative of v<i in 
Equation (3.6) with respect to time. The result is 



dVd dv d dd vr e [r e r s sin 2 0 - (r% + r 5 2 )sin0 + r e r s ] dd 
dt dO dt (rj + r] - 2r e r s sin 0) 3 / 2 dt 

In deriving this equation, the result of Equation (3.7) is used. The result of 
this equation is shown in Figure 3.5 and the maximum rate of change of the 
frequency occurs at 0 = 7r/2. 

The corresponding maximum rate of change of the speed is 



dVd vr e dO/dt 

dt max ” yjr l + rj - 2 r e r x 



~ 0.178 m/s 2 



0 7r/2 



(3.15) 
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FIGURE 3.5 The rate of change of the speed versus angle 6 . 



In this equation, only the magnitude is of interest, thus, the sign is neglected. 
The corresponding rate of change of the Doppler frequency is 



. dv d f r 0.178 x 1575.42 x 10 6 
8fdA ™ = — = 3^08 



- 0.936 Hz/s 



This value is also very small. If the frequency accuracy measured through the 
tracking program is assumed on the order of 1 Hz, the update rate is almost 
one second, even at the maximum Doppler frequency changing rate. 



3.9 RATE OF CHANGE OF THE DOPPLER FREQUENCY DUE TO USER 
ACCELERATION 

From the previous two sections, it is obvious that the rate of change of the 
Doppler frequency caused by the satellite motion is rather low; therefore, it 
does not affect the update rate of the tracking program significantly. 

Now let us consider the motion of the user. If the user has an acceleration 
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of 1 g (gravitational acceleration with a value of 9.8 m/s 2 ) toward a satellite, 
the corresponding rate of change of the Doppler frequency can be found from 
Equation (3.15) by replacing dvjdt by g. The corresponding result obtained 
from Equation (3.16) is about 51.5 Hz/s. For a high-performance aircraft, the 
acceleration can achieve several g values, such as 7 g. The corresponding rate of 
change of the Doppler frequency will be close to 360 Hz/s. Comparing with the 
rate of change of the Doppler frequency caused by the motions of the satellite 
and the receiver, the acceleration of the receiver is the dominant factor. 

In tracking the GPS signal in a software GPS receiver two factors are used 
to update the tracking loop: the change of the carrier frequency and the align- 
ment of the input and the locally generated C/A codes. As discussed in Section 
3.5, the input data adjustment rate is about 20 ms due to the Doppler frequency 
on the C/A code. If the carrier frequency of the tracking loop has a band- 
width of the order of 1 Hz and the receiver accelerates at 7 g, the tracking 
loop must be updated approximately every 2.8 ms (1/360) due to the carrier 
frequency change. This might be a difficult problem because of the noise in 
the received signal. The operation and performance of a receiver tracking loop 
greatly depends on the acceleration of the receiver. 



3.10 KEPLER'S LAWS* 11 ' 12 ) 

In the previous section, the position of a satellite is briefly described. This infor- 
mation can be used to determine the differential power level and the Doppler 
frequency on the input signal. However, this information is not sufficient to 
calculate the position of a satellite. To find the position of a satellite, Kepler’s 
laws are needed. The discussion in this section provides the basic equations to 
determine a satellite position. 

Kepler’s three laws are listed below (see Chapter 1 in ref. 11): 

First Law: The orbit of each planet is an ellipse with the sun at a focus. 

Second Law: The line joining the planet to the sun sweeps out equal areas 
in equal times. 

Third Law: The square of the period of a planet is proportional to the cube 
of its mean distance from the sun. 

These laws also apply to the motion of the GPS satellites. The satellite orbit 
is elliptical with the earth at one of the foci. Figure 3.6 shows the orbit of a GPS 
satellite. The center of the earth is at F and the position of the satellite is at S. 
The angle v is called the actual anomaly. In order to illustrate the basic concept, 
the shape of the ellipse is overemphasized. The actual orbit of the satellite is 
very close to a circle. The point nearest to the prime focus is called the perigee 
and the farthest point is called the apogee. 

Kepler’s second law can be expressed mathematically as (Figure 3.6) 
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t - t p 
Ai 



T 

ira s b s 



(3.17) 



where t presents the satellite position at time t, t p is the time when the satellite 
passes the perigee, Ai is the area enclosed by the lines t = t, t = t p , and the 
ellipse, T is the period of the satellite, a s and b s are the semi-major and semi- 
minor axes of the orbit, and it a s b s is the total area of the ellipse. This equation 
states that the time to sweep the area A } is proportional to the time T to sweep 
the entire area of the ellipse. 

The third law can be stated mathematically as 




where n = GM = 3.986005 x 10 14 meters 3 /sec 2 (ref. 12) is the gravitational 
constant of the earth. Thus, the right-hand side of this equation is a constant. In 
this equation the semi-major axis a s is used rather than the mean distance from 
the satellite to the center of the earth. In reference 11 it is stated that a s can be 
used to replace the mean distance because the ratio of a s to the mean distance 
r s is a constant. This relationship can be shown as follows. If one considers the 
area of the ellipse orbit equal to the area of a circular orbit with radius r s , then 
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ira s b s - 7rr 2 or — = (3.19) 

r s b s 

Since a s , b s , r s are constants, a s and r s is related by a constant. 



3.11 KEPLER'S EQUATION! 11 ' 13 ) 

In the following paragraphs Kepler’s equation will be derived and the mean 
anomaly will be defined. The reason for this discussion is that the information 
given by the GPS system is the mean anomaly rather than the actual anomaly 
that is used to calculate the position of a satellite. 

In order to perform this derivation, a few equations from the previous chapter 
will be repeated here. The eccentricity is defined as 



where c s is the distance from the center of the ellipse to a focus. For an ellipse, 
the e s value is 0 < e s < 1. When a s —b s , then e s — 0, which represents a circle. 
The eccentricity e s can be obtained from data transmitted by the satellite. 

In Figure 3.7 an elliptical satellite orbit and a fictitious circular orbit are 
shown. The center of the earth is at F and the satellite is at S. The area A\ is 
swept by the satellite from the perigee point to the position S. This area can be 
written as 



A\ = area PSV- A 2 (3.21) 

In the previous chapter Equation (2.24) shows that the heights of the ellipse 
and the circle can be related as 



QP _ a s 
~SP ~b 7 



(3.22) 



Therefore, the area PSV can be obtained from area PQV as 



= — — a~E a 1 sin E cos E = - 



—(E - sin E cos E) 



(3.23) 




FIGURE 3.7 Fictitious and actual orbits. 



where the angle E is called eccentric anomaly. The area of triangle A 2 is 

A 2 = — SPxPF=— — QP X PF = — —a s sin E(c s - a s cos E ) 

2 2 a s 2 a s 

= sin E(e s a s - a s cos E ) = Us ^ s ( e s sin E - sin £ cos E) (3.24) 

In the above equation, the relation in Equation (3.20) is used. Substituting 
Equations (3.23) and (3.24) into (3.21) the area is 

Ai - a ^- ( E - e s sin E) (3.25) 

Substituting this result into Equation (3.17), Kepler’s second law, the result is 
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The next step is to define the mean anomaly M and from Equation (3.26) 
the result is 



M = (E 



sin E) 




(3.27) 



If one defines the mean motion n as the average angular velocity of the satellite, 
then from Equation (3.18) the result is 




Substituting this result into Equation (3.27) the result is 

M = (E - e s sin E) = n(t - t p ) (3.29) 

This is referred to as Kepler’s equation. From this equation one can see that M 
is linearly related to f, therefore, it is called the mean anomaly. 



3.12 TRUE AND MEAN ANOMALY 

The information obtained from a GPS satellite is the mean anomaly M. From 
this value, the true anomaly must be obtained because the true anomaly is used 
to find the position of the satellite. The first step is to obtain the eccentric 
anomaly E from the mean anomaly, Equation (3.29) relates M and E. Although 
this equation appears very simple, it is a nonlinear one; therefore, it is difficult 
to solve analytically. This equation can be rewritten as follows: 

E = M + e s sin E (3.30) 

In this equation, e s is a given value representing the eccentricity of the satel- 
lite orbit. Both e s and M can be obtained from the navigation data of the satel- 
lite. The only unknown is E. One way to solve for E is to use the iteration 
method. A new E value can be obtained from a previous one. The above equa- 
tion can be written in an iteration format as 

E i+l =M + e s smEi (3.31) 

where E i+1 is the present value and E, is the previous value. One common 
choice of the initial value of E is E 0 = M. This equation converges rapidly 
because the orbit is very close to a circle. Either one can define an error signal 
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as E err = E i+l - E t and end the iteration when E err is less than a predetermined 
value, or one can just iterate Equation (3.31) a fixed number of times (i.e., from 
5 to 10). 

Once the E is found, the next step is to find the true anomaly v. This relation 
can be found by referring to Figure 3.7. 



OP c s - PF 



(3.32) 



Now let us find the distance r in terms of angle v. From Figure 3.6, applying 
the law of cosine to the triangle GSF, the following result is obtained 

r' 1 r 2 + 4rc s cos v + 4c s 2 (3.33) 

where r and r are the distance from the foci G and F to the point S. For an 
ellipse, 



r + r = 2 a s 

Substituting this relation into Equation (3.33), the result is 

al-c] o f (l - e 2 s ) 
a s + c s cos v 1 + e s cos v 

Substituting this value of r into Equation (3.32) the result is 



(3.34) 



cos E = 



e s + cos v 
1 + e s cos v 



Solve for v and the result is 



cos E - e s 
1 - e s cos E 



(3.36) 



(3.37) 



This solution generates multiple solutions for v because cos v is a multival- 
ued function. One way to find the correct value of v is to keep these angles E 
and v in the same half plane. From Figure 3.7 one can see that the angles E 
and v are always in the same half plane. 

Another approach to determine v is to find the sin v. (i3 > If one takes the 
square on both sides of the above equation, the result is 
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Thus, the true anomaly can be found as 



v V] sign (v 2 ) 



(3.41) 



where sign(j» 2 ) provides the sign of P 2 \ therefore, it is either +1 or -1. It is 
interesting to note that to find the true anomaly only M and e s are needed. 
Although the semi-major axis a s appears in the derivation, it does not appear 
in the final equation. 



3.13 SIGNAL STRENGTH AT USER LOCATIONfbS^-ie) 

In this section the signal strength at the user location will be estimated. The 
signal strength can be obtained from the power of the transmitting antenna, the 
beam width of the antenna, the distance from the satellite to the user, and the 
effective area of the receiving antenna. The power amplifier of the transmitter 
is 50 w® (or 17 dBw). The input to the transmitting antenna is 14.3 dBw. (8 J 
The difference might be due to impedance mismatch or circuit loss. 

The gain of the transmitting antenna can be estimated from the beam width 
(or solid angle) of the antenna. The solid angle is denoted as 9, which is 21.3 
degrees. The area on the surface of a sphere covered by the angle 9 can be 
obtained from Figure 3.9 as 




FIGURE 3.9 Area facing solid angle 6 . 
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2-k(t sin 6)rd0 = 2 -k r 



- 2irr 2 (- cos 0 )|q - 27rr 2 (l - cos 6 



The ratio of this area to the area of the sphere can be considered as the gain of 
the transmitting antenna, which can be written as 



4tt r 2 



27rr 2 (l - cos ^)| 21 . 3 ° 0.683 



= 29.28 = 14.7 dB 



Using 14.3 dBw as the input to the antenna, the output of the antenna should 
be 29 dBw (14.3 + 14.7). However, the transmitting power level is listed as 
478.63 w,0 4 45) which corresponds to 26.8 dBw. This difference between the 
power levels might be due to efficiency of the antenna and the accuracy of the 
solid angle of the antenna because the power cannot be cut off sharply at a 
desired angle. 

If the receiving antenna has a unit gain, the effective area is (16 ) 



A 



(3.44) 



where X is the wavelength of the receiving signal. 

The received power is equal to the power density multiplied by the effective 
area of the receiving antenna. The power density is equal to the radiating power 
divided by the surface of the sphere. The receiving power can be written as 

47 vR~ u 47 rR 2 , 47T (47 rR su ) 2 



where R su is the distance from the satellite to the user. Assume R su = 25785 x 
10 3 m, which is the farthest distance as shown in Figure 3.1. Using 478.63 W 
as the transmitting antenna and the wavelength X = 0.19 m, the receiving power 
P r calculated from the above equation is 1.65 x 10 16 w (or -157.8 dBw). If 
the loss through the atmosphere is taken into consideration, the received power 
is close to the minimum required value of -160 dBw. 

The power level at the receiver is shown in Figure 3.10. It is a function 
of the elevation angle/') At zenith and horizon, the powers are at -160 dBw. 
The maximum power level is -158 dBw, which occurs at about 40 degrees. If 
the receiving antenna is taken into consideration, the received power will be 
modified by its antenna pattern. 
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FIGURE 3.10 Power level versus elevation angle. 



3.14 SUMMARY 

This chapter discusses the orbits of the GPS satellite. The orbit is elliptical but 
it is very close to a circle. Thus, the circular orbit is used to figure the power 
difference to the receiver and the Doppler frequency shift. This information is 
important for tracking the satellite. In order to find the position of a satellite the 
actual elliptical satellite orbit must be used. To discuss the motion of the satellite 
in the elliptical-shaped orbit, Kepler’s laws are introduced. Three anomalies are 
defined: the mean M, the eccentric E, and the true v anomalies. Mean anomaly 
M and eccentricity e s are given from the navigation data of the satellite. Eccen- 
tric anomaly E can be obtained from Equation (3.30). True anomaly v can be 
found from Equations (3.40) and (3.41). Finally, the receiving power at the user 
location is estimated. 
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CHAPTER. FOUR 



Earth-Centered, Earth-Fixed 
Coordinate System 



4.1 INTRODUCTION 

In the previous chapter the motion of the satellite is briefly discussed. The true 
anomaly is obtained from the mean anomaly, which is transmitted in the navi- 
gation data of the satellite. In all discussions, the center of the earth is used as 
a reference. In order to find a user position on the surface of the earth, these 
data must be related to a certain point on or above the surface of the earth. 
The earth is constantly rotating. In order to reference the satellite position to a 
certain point on or above the surface of the earth, the rotation of the earth must 
be taken into consideration. This is the goal of this chapter. 

The basic approach is to introduce a scheme to transform the coordinate sys- 
tems. Through coordinate system transform, the reference point can be moved 
to the desired coordinate system. First the direction cosine matrix, which is used 
to transform from one coordinate system to a different one, will be introduced. 
Then various coordinate systems will be introduced. The final transform will 
put the satellite in the earth-centered, earth-fixed (ECEF) system. Finally, some 
perturbations will be discussed. The major portion of this discussion is based 
on references 1 and 2. 

In order to perform the transforms, besides the eccentricity e s and mean 
anomaly M, additional data are obtained from the satellite. They are the semi- 
major of the orbit a s , the right ascension angle fi , the inclination angle i, and 
the argument of the perigee co. Their definitions will also be presented in this 
chapter. 
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4.2 DIRECTION COSINE MATRIX* 1 " 3 ) 

In this section, the direction cosine matrix will be introduced. A simple two- 
dimensional example will be used to illustrate the idea, which will be extended 
into a three-dimensional one without further proof. Figure 4. 1 shows two two- 
dimensional systems (xi, yi) and (X2, y 2 )- The second coordinate system is 
obtained from rotating the first system by a positive angle a. A point p is used 
to find the relation between the two systems. The point p is located at (X \ , 
Fi) in the (x \ , yi) system and at (X 2 , Y 2 ) in the (X 2 , y 2 ) system. The relation 
between ( X 2 , Y 2 ) and {X u Y{) can be found from the following equations: 



X 2 = X\ cos a + Y\ sin a = Xi cos(Xi on X 2 ) + Y\ cos(Ti on X 2 ) 

Y 2 - -Xi sin a + Y x cos a = Xi cos(Xi on Y 2 ) + Y x cos(Ti on Y 2 ) (4.1) 



In matrix form this equation can be written as 




FIGURE 4.1 Two coordinate systems. 
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' X 2 1 [ cos(Xi on X 2 ) cos(yi on X 2 ) Xj 

Y 2 “ cos(Xi on Y 2 ) cos(Fi on Y 2 ) Ti 



The direction cosine matrix is defined as 



cos(Xi onX 2 ) cos(TionX 2 ) 
cos(Xi on Y 2 ) cos(Fi on Y 2 ) 



This represents that the coordinate system is transferred from system 1 to sys- 
tem 2. 

In a three-dimensional system, the directional cosine can be written as 



cl 



cos(Xi on X 2 ) 
cos^i on T 2 ) 
cos(Xi on Z 2 ) 



cos(yi on X 2 ) 
cos(yi on y 2 ) 
cos(yi on Z 2 ) 



cos(Zi on X 2 ) ' 
cos(Zi on y 2 ) 
cos(Zi on Z 2 ) . 



(4.3) 



Sometimes it is difficult to make one single transform from one coordinate 
to another one, but the transform can be achieved in a step-by-step manner. For 
example, if the transform is to rotate angle a. around the /-axis and rotate angle 
1 6 around the y-axis, it is easier to perform the transform in two steps. In other 
words, the directional cosine matrix can be used in a cascading manner. The 
first step is to rotate a positive angle a around the z-axis. The corresponding 
direction cosine matrix is 



I cos a sin a 0 1 

-sin a cos a 0 (4.4) 

0 0 lj 

The second step is to rotate a positive angle j8 around the x-axis; the corre- 
sponding direction cosine matrix is 



C\ - 0 cos /3 sin j8 (4.5) 

.0 -sin j8 cos/S. 



The overall transform can be written as 
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1 0 


0 ■ 


' cos a 


sin a 0 ] 


0 cos |8 


sin |3 


-sin a. 


cos a 0 


0 -sin j3 


cos /3. 


. 0 


0 lJ 


cos a 


sin a 


° 1 


-sin a cos j£ 


’ cos a cos j8 


sin j3 


sin a sin (3 


-cos 


a sin |3 


cos (3J 



It should be noted that the order of multiplication is very important; if the order 
is reversed, the wrong result will be obtained. 

Suppose one wants to transform from coordinate system 1 to system n 
through system 2, 3, ... n - 1 . The following relation can be used: 



C n l =C n „_ 1 ---C\cl (4.7) 

In general, each C\ , represents only one single transform. This cascade method 
will be used to obtain the earth-centered, earth-fixed system. 



4.3 SATELLITE ORBIT FRAME TO EQUATOR FRAME TRANSFORM^) 

The coordinate system used to describe a satellite in the previous chapter can 
be considered as the satellite orbit frame because the center of the earth and 
the satellite are all in the same orbit plane. Figure 4.2 shows such a frame, and 
the x-axis is along the direction of the perigee and the z-axis is perpendicular 
to the orbit plane. The y-axis is perpendicular to the x and z axes to form a 
right-hand coordinate system. The distance r from the satellite to the center of 
the earth can be obtained from Equation (3.35) as 



q,(l - # 

1 + e s cos v 



(4.8) 



where a s is the semi-major of the satellite orbit, e s is the eccentricity of the satel- 
lite orbit, v is the true anomaly, which can be obtained from previous chapter. 
The value of cos v can be obtained from Equation (3.37) as 



(4.9) 



where E is the eccentric anomaly, which can be obtained from Equation (3.30). 

Substituting Equation (4.9) into Equation (4.8) the result can be simplified 
as 
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r — 0,(1 - e s cos E) (4.10) 

The position of the satellite can be found as 

x — r cos v 
y = r sin y 

z = 0 (4.11) 

This equation does not reference any point on the surface of the earth but refer- 
ences the center of the earth. It is desirable to reference to a user position that 
is a point on or above the surface of the earth. 

First a common point must be selected and this point must be on the surface 
of the earth as well as on the satellite orbit. The satellite orbit plane intercepts 
the earth equator plane to form a line. An ascending node is defined along 
this line toward the point where the satellite crosses the equator in the north 
(ascending) direction. The angle go between the perigee and ascending node in 
the orbit plane is referred to as the argument of the perigee. This angle infor- 
mation can be obtained from the received satellite signal. Now let us change 
the x-axis from the perigee direction to the ascending node. This transform can 
be accomplished by keeping the z-axis unchanged and rotating the x-axis by 
the angle go as shown in Figure 4.3. In Figure 4.3 the y-axis is not shown. The 
x,-axis and the z,-axis are perpendicular and the y,-axis is perpendicular to the 
XfZi plane. The corresponding direction cosine matrix is 
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equinox 

FIGURE 4.3 Earth equator and orbit plane. 



I cos co -sin co 0] 

sin co cos co 0 (4.12) 

0 0 lJ 

In this equation the angle co is in the negative direction; therefore the sin co 
has a different sign from Equation (4.4). This rotation changes the xi-axis to 
^2-axis. 

The next step is to change from the orbit plane to the equator plane. This 
transform can be accomplished by using the X 2 -axis as a pivot and rotate angle 
/. This angle i is the angle between the satellite orbit plane and the equator 
plane and is referred to as the inclination angle. This inclination angle is in the 
data transmitted by the satellite. The corresponding direction cosine matrix is 



C\= 0 cos i -sin/ (4.13) 

L 0 sin i cos i J 

The angle i is also in the negative direction. After this transform, the Z 3 -axis is 
perpendicular to the equator plane rather than the orbit of the satellite and the 
x 3 -axis is along the ascending point. 

There are six different orbits for the GPS satellites; therefore, there are six 
ascending points. It is desirable to use one x-axis to calculate all the satellite 
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positions instead of six. Thus, it is necessary to select one x-axis; this subject 
will be discussed in the next section. 



4.4 VERNAL EQUINOX< 2 ) 

The vernal equinox is often used as an axis in astrophysics. The direction of the 
vernal equinox is determined by the orbit plane of the earth around the sun (not 
the satellite) and the equator plane. The line of intersection of the two planes, 
the ecliptic plane (the plane of the earth’s orbit) and the equator, is the direction 
of the vernal equinox as shown in Figure 4.4. 

On the first day of spring a line joining from the center of the sun to the 
center of the earth points in the negative direction of the vernal equinox. On 
the first day of autumn a line joining from the center of the sun to the center 
of the earth points in the positive direction of the vernal equinox as shown in 
Figure 4.5. 

The earth wobbles slightly and its axis of rotation shifts in direction slowly 
over the centuries. This effect is known as precession and causes the line-of- 
intersection of the earth’s equator and the ecliptic plane to shift slowly. The 
period of the precession is about 26,000 years, so the equinox direction shifts 
westward about 50 (360 x 60 x 60/26000) arc-seconds per year and this is a 
very small value. Therefore, the vernal equinox can be considered as a fixed 
axis in space. 

Again referring to Figure 4.3, the A' 3 -axis of the last frame discussed in the 
previous section will be rotated to the vernal equinox. This transform can be 
accomplished by rotating around the Z 3 -axis an angle referred to as the right 
ascension. This angle is in plane of the equator. The direction cosine matrix is 
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I cos U -sin Q 0] 

sin fi cos 0 0 (4. 14) 

0 0 lj 

This last frame is often referred to as the earth-centered inertia (ECI) frame. 
The origin of the ECI frame is at the earth’s center of mass. In this frame the 
Z 4 -axis is perpendicular to the equator and the jo-axis is the vernal equinox 
and in the equator plane. This frame does not rotate with the earth but is fixed 
with respect to stars. In order to reference a certain point on the surface of the 
earth, the rotation of the earth must be taken into consideration. This system is 
referred to as the earth-centered, earth-fixed (ECEF) frame. 



4.5 EARTH ROTATION^) 

In this section two goals will be accomplished. The first one is to take care 
of the rotation of the earth. The second one is to use GPS time for the time 
reference. 

First let us consider the earth rotation. Let the earth turning rate be fl ie and 
define a time t er such that at t er = 0 the Greenwich meridian aligns with the 
vernal equinox. The vernal equinox is fixed by the Greenwich meridian rotates. 
Referring to Figure 4.6, the following equation can be obtained 




62 



EARTH-CENTERED, EARTH-FIXED COORDINATE SYSTEM 




Meridian at 
timet^ 

FIGURE 4.6 Rotation of the earth. 



Q er -U - (l le t er (4.15) 

where 0 er is the angle between the ascending node and the Greenwich meridian, 
the earth rotation rate Q ie = 7.2921151467 x 10 5 rad/sec. When t er = 0, Q er = 
1] , this means that the Greenwich meridian and the vernal equinox are aligned. 

If the angle U er is used in Equation (4.14) to replace 0, the x-axis will 
be rotating in the equator plane. This x-axis is the direction of the Greenwich 
meridian. Using this new angle in Equation (4.14) the result is 

I COS tier -sin Q, er 0] 

sinQ^r cos Q er 0 (4.16) 

0 0 lJ 

In this equation the rotation of the earth is included, because time is included 
in Equation (4.15). Using this time t er in the system, every time the Greenwich 
meridian is aligned with the vernal equinox, t er = 0. The maximum length of this 
time is a sidereal day, because the Greenwich meridian and the vernal equinox 
are aligned once every sidereal day. 

The time t er should be changed into the GPS time t. The GPS time t starts 
at Saturday night at midnight Greenwich time. Thus, the maximum GPS time 
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is seven solar days. It is obvious that the time base t er and the GPS time t are 
different. A simple way to change the time t er to GPS time t is a linear shift 
of the time base as 



(4.17) 



where At can be considered as the time difference between the time based on 
t er and the GPS time t. Substituting this equation into Equation (4.15), the result 
is 



Q er = U - Q ie t er = Q - Q ie t - Q ie A t = 
where fl e = fi - a and a = 0 ie A t 



-Q ie t = Q e - Q ie t 

(4.18) 



The reason for changing to this notation is that the angle 0 - a is considered as 
one angle fi e , and this information is given in the GPS ephemeris data. How- 
ever, this relation will be modified again in Section 4.7 and the final result will 
be used to find U er in Equation (4.16). Before the modification of Q e , let us 
first find the overall transform. 



4.6 OVERALL TRANSFORM FROM ORBIT FRAME TO EARTH-CENTERED, 
EARTH-FIXED FRAME 

In order to transform the positions of the satellites from the satellite orbit frame 
to the ECEF frame, there need to be two intermediate transforms. The overall 
transform can be obtained from Equation (4.7). Substituting the results from 
Equations (4.16), (4.13), and (4.12) into (4.7), the following result is obtained: 



= C^C^Cj 2 



COS Q. er 
sin Q er 
0 




0 

0 



COS Qer 
sin Q er 
0 



-sin Q er cos i 
cos Qer COS i 



sin Q er sin i 
-cos Q er sin i 



0 




64 EARTH-CENTERED, EARTH-FIXED COORDINATE SYSTEM 



cos Q e r cos co - sin Q er cos i sin 01 
sin Q er cos 01 + cos Q er cos i sin a 



i - sin fi er 



sin Q er sin i 
-cos Q er sin 1 






r cos Q er cos(v + 01) - r sin Q er cos i sin(V + 01) ] 

r sin Q er cos (y + oi) + r cos Q er cos i sinO + 01) ( 4 . 19 ) 

rsinisin( V +o>) \ 

This equation gives the satellite position in the earth-centered, earth-fixed coor- 
dinate system. 

In order to calculate the results in the above equation, the following data are 
needed: (1) a s : semi-major axis of the satellite orbit; (2) M: mean anomaly; (3) 
e s : eccentricity of the satellite orbit; (4) i: inclination angle; (5) co: argument of 
the perigee; (6) 0 -a: modified right ascension angle; (7) GPS time. The first three 
constants are used to calculate the distance r from the satellite to the center of the 
earth and the true anomaly v as discussed in Section 3.12. The three values i, co, 
and fi -a are used to transform from the satellite orbit frame to the ECEF frame. 
In order to find U er in the above equation the GPS time is needed. 



4.7 PERTURBATIONS 

The earth is not a perfect sphere and this phenomenon affects the satellite orbit. 
In addition to the shape of the earth, the sun and moon also have an effect on 
the satellite motion. Because of these factors the orbit of the satellite must be 
modified by some constants. The satellites transmit these constants and they 
can be obtained from the ephemeris data. 

Equation (4.19) is derived based on the assumption that the orbit of the satel- 
lite is elliptical; however, the orbit is not a perfect elliptic. Thus, the parameters 
in the equations need to be modified. This section presents the results of the 
correction terms. 

In Equation (4.15) the right ascension Q will be modified as 

(t-t oe ) (4.20) 

where t is the GPS time, t oe is the reference time for the ephemeris, and 0 is 
the rate of change of the right ascension. In this equation it is implied that the 
right ascension is not a constant, but changes with time. The ephemeris data 
transmitted by the satellite contain t oe and 0. Substituting this equation into 
Equation (4.18), the result is 

Q er = Q -a + Q(t- t oe ) - 12 iet = II e + 0 (t ~ t oe ) - 0 ie t 



(4.21) 
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where 0 e is contained in the ephemeris data. 

The mean motion in Equation (3.28) must be modified as 

nm*i + An = J-^ + An (4.22) 

V a s 

where An is the correction term that is contained in the ephemeris data. The 
mean anomaly must be modified as 

M - Mo + n(t - t oe ) (4.23) 

where Mq is the mean anomaly at reference time, which can be obtained from 

the ephemeris data. This value M will be used to find the true anomaly v. 

There are six constants C us , C uc , C rs , C rc , C- ls , and C ic and they are used to 
modify v + co, r, and i in Equation (4.19) respectively. Let us introduce a new 
variable 0 as 

<t>=v + co (4.24) 

The correction term to v + co is 

6(v + o) = kt> = C us sin 2<t> + C uc cos 20 (4.25) 

and the new v + co is 

v + co => v + co + 5(v + o}) (4.26) 

The correction to distance r is 

dr = C rs sin 20 + C rc cos 20 (4.27) 

and the new r is 

r^r + br (4.28) 

The correction to inclination i is 

bi - C is sin 20 + C, c cos 20 (4.29) 

and the new inclination i is 

i?4'i + bi (4.30) 
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Substituting these new values into Equation (4.19) will produce the desired 
results. 



4.8 CORRECTION OF GPS SYSTEM TIME AT TIME OF 
TRANSMISSION^' 6 ) 

In Equations (4.21) and (4.23) the GPS time is used, and this time is often referred 
to at the time of transmission. (This section discusses only the correction of this 
time. The actual obtaining of the time of transmission will be discussed in Sec- 
tion 9.10.) The signals from the satellites are transmitted at the same time except 
for the clock error in each satellite. The time of receiving t u is the time the signal 
arrives at the receiver. The relation between the t and t u is 



t u - t + Pi/C 

t = t„- Pt/c (4.31) 



where p, is the pseudorange from satellite i to the receiver and c is the speed of 
light. Since the pseudorange from each satellite to the receiver is different, the 
time of receiving is different. However, in calculating user position, one often 
uses one value for time. The time of receiving t u is a reasonable selection. If 
a time of receiving t u is used as a reference, from the above equation the time 
of transmission from various satellites is different. The time of transmission is 
the receiving time minus the transit time. This time is represented by t and is 
referred to as the time of transmission corrected for the transit time. 

The t value must be corrected again from many other factors. However, in 
order to correct t, the t value must first be known. This requirement makes the 
correction process difficult. To simplify this process, let t c represent the coarse 
GPS system time at time of transmission corrected by transit time. The value 
t c can be obtained from time of the week (TOW), which will be presented in 
Section 5.9. For the present discussion, let us assume that the t c value is already 
obtained. The time t k shall be the actual total time difference between the time 
t c and the epoch time t oe and must account for the beginning or end of the week 
crossovers. That is, the following adjustments must be made on t c : 

If t k - t c - t oe > 302400 then t k = t k - 604800 or t c => t c - 604800 
If t k - t c - t oe < -302400 then t k =t k + 604800 or t e s* t c + 604800 

(4.32) 

where t oe can be obtained from ephemeris data, 302,400 is the time of half a 
week in seconds. The time of a week in seconds is 604,800 (7 x 24 x 3600). 

The following steps can be used to correct the GPS time t. From Equation 
(4.22), the mean motion is calculated as 
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(4.33) 



where g, = 3.986005 x 10 14 meters 3 /sec 2 is the earth’s universal gravitational 
parameter and is a constant, and An are obtained from ephemeris data. 
From this n value the mean anomaly can be found from Equation (4.23) as 

M = M 0 + n(t c - t oe ) (4.34) 



where Mq is in the ephemeris data. In this equation t c is used instead of t as t 
is not derived yet. 

The eccentric anomaly E can be found from Equations (3.29) or (3.30) 
through iteration as 



E = M + e s sin E (4.35) 

where e s is eccentricity of the satellite orbit, which can be obtained from the 
ephemeris data. Let us define a constant F as 



F = — = -4.442807633 x 10 10 sec/(meter) 1 / 2 (4.36) 
c l 

where /r is the earth’s universal gravitational parameter and c is the speed of 
light. The relativistic correction term is 



A t r = Fe s ■sfa s sin E (4.37) 

The overall time correction term is 



At - afo + a.fi(t c - t oc ) + af2(t c - t oc ) 2 + A t r -T GD (4.38) 

where T GD , t oc , afo, aj\ , aj 2 are clock correction terms and T GD is to account 
for the effect of satellite group delay differential. They can be obtained in the 
ephemeris data. The GPS time of transmission t corrected for transit time can 
be corrected as 



(4.39) 



This is the time t that will be used for the following calculations. 




68 EARTH-CENTERED, EARTH-FIXED COORDINATE SYSTEM 



4.9 CALCULATION OF SATELLITE POSITION* 5 ' 6 ) 

This section uses all the information from the ephemeris data to obtain a satellite 
position in the earth-centered, earth-fixed system. These calculations require the 
information obtained from both Chapters 3 and 4; therefore, this section can be 
considered as a summary of the two chapters. 

Equation (4.19) is required to calculate the position of the satellite. In this 
equation there are five known quantities: r, v + o>, i, and Q er . These quantities 
appear on the right side of the equation and the results represent the satellite 
position. Let us find these five quantities. 

First let us find the value of r from Equation (4.10) as 

r = 0,(1 - e s cos E) (4.40) 

In this equation, the value E must be calculated first from ephemeris data. In 
order to find the r value the following steps must be taken: 

1. Use Equation (4.22) to calculate n where ji is a constant; a s and An can 
be obtained from the ephemeris data. 

2. Use Equation (4.34) to calculate M where M 0 and t oe can be obtained from 
ephemeris data and t c can be obtained from the discussion in Section 9.10. 

3. The value of E can be found from Equation (4.35), where e s can be 
obtained from the ephemeris data. The iteration method will be used in 
this operation. 

4. Once E is obtained, the value of r can be found from Equation (4.40). 

In the above four steps, the first three steps are to find the value of E. Once 
E is calculated, Equations (4.36)-(4.39) can be used to find the corrected GPS 
time t. 

Now let us find the true anomaly v. This value can be found from Equations 
(3.40) and (3.41) as 



V 2 = sin 
v = v\ sign(i> 2 ) 



1 cos E- e s \ 
v 1 - e?cos E ) 



/l - e 2 s sin E 
1 - e s cos E 



(4.41) 



The argument w can be found from the ephemeris data. Using the definition 
in Equation (4.24), the value of 0 is 
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The following correction terms are needed 

8<j) - C us sin 2 4> + C uc cos 20 
dr - C rs sin 2 <f) + C rc cos 20 

8i = Ci S sin 2<j> + Ci c cos 20 (4.43) 

where the C us , C uc , C rs , C rc , Ci S , Ci r are from ephemeris data: 

0 => 0 + 50 

r => r + <5r (4.44) 

The inclination angle i can be obtained from the ephemeris data and be corrected 
as 

i => i + 8i + idot(t - t oe ) (4.45) 

where idot can be obtained from the ephemeris data. The last term to be found 
is 



Q er = Q e + to(t-t oe )-Q ie t (4.46) 

where the earth rotation rate Q is a constant, 0 , and t oe are obtained from 

the ephemeris data. It should be noted that the corrected GPS time t is used in 
the above two equations. 

Once all the necessary parameters are obtained, the position of the satellite 
can be found from Equation (4.19) as 



r cos Q er cos cj> - r sin fl er cos i sin <j) ' 
r sin er cos 4>+r cos Q er cos i sin <£ 
r sin i sin 0 



(4-47) 



The satellite position calculated in this equation is in the ECEF frame. There- 
fore the satellite position is a function of time. From the x, y, z, and the pseudo- 
range of more than four satellites the user’s position can be found from results 
in Chapter 2. The actual calculation of the pseudorange is discussed in Sec- 
tion 9.9. 



4.10 COORDINATE ADJUSTMENT FOR SATELLITES 

Using the earth-centered, earth-fixed coordinate system implies that the earth’s 
rotation is taken into consideration. The satellite position calculated from Sec- 
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tion 4.9 is based on the GPS time of transmission t corrected for transit time. 
However, the user position is calculated at the time of receiving. Since the satel- 
lite and user positions are calculated at different times, they are in different 
coordinate systems. This will cause an error in the user position. As discussed 
in Section 3.3., if the user is on the equator of the earth and an approximate 
signal traveling time of 76 ms is assumed, the user position is moved about 36 
m (2ir x 6368 x 10 3 x 76/(24 x 3600 x 10 3 )) due to the rotation of the earth. 

In order to obtain the correct user position, a single coordinate system should 
be used. Since the user position is measured at the time of receiving, it is appro- 
priate to use this time in the coordinate system. The satellite position calculated 
should be referenced to this time. This correction does not mean to change the 
satellite position, but only changes the coordinate system of the satellites. 

In order to reference the GPS time at the time of receiving, the coordinate 
system of each satellite must be separately modified. Using the time of receiving 
as reference, the time of transmission of satellite i in the new coordinate system 
is the time of receiving minus the transition time as shown in Equation (4.31). 
The transit time cannot be determined before the user position is calculated 
because of the unknown user clock bias. Only when the user position is obtained 
can the pseudorange be found. Once the user position is found, the pseudorange 
can be found from the user position to the satellite position. The earth rotation 
appears only in Equation (4.46). Equations (4.46) and (4.47) are used in the 
operation. 

The following steps can be taken to improve user position accuracy. 

1. From the satellite and user position the transit time t t can be found as 



h = V( x-x u ) 2 + (y-y u ) 2 + (z-z u ) 2 /c (4.48) 

where x, y, z, and x u , y u , z u are the coordinates of the satellite and the 
user respectively, c is the speed of light. 

2. Use the transit time to modify the angle U er in Equation (4.46) as 

£l er ^Q er -{l ie t t (4.49) 

3. Use the new value of fi er in Equation (4.47) to calculate the position of 
the satellite x, y, and z. 

4. The above operations should be performed on every satellite. From these 
values a new user position x u ,y u , z u will be calculated. 

5. Repeat steps 1, 2, 3, and 4 again to obtain a new set of x, y, and z. When 
the old and new sets are within a predetermined value, the new set can be 
considered as the position of the satellite in the new coordinate system. 
It usually requires calculating the x, y, and z values only twice. 

6. These new x, y, and z values will be used to find the user position. 
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4.11 EPHEMERIS DATA( 4 - 6 ) 

In the previous sections several constants and many ephemeris data are used 
in the calculations. This section lists all these constants and the ephemeris data 
used in the calculations. The details of the ephemeris data transmitted by the 
satellites will be presented in the next chapter. 

The constants are listed as follows:' 4 * 

/a = GM = 3.986005 x 10 14 meters 3 /sec 2 , which is the WGS-84 value of the 
earth’s universal gravitational parameter. 

£l ie - 7.2921151467 x 10 5 rad/sec, which is the WGS-84 value of the 
earth’s rotational rate. 

7T = 3.1415926535898. 

c = 2.99792458 x 10 8 meter/sec, which is the speed of light. 

The ephemeris data are: 

M 0 : mean anomaly at reference time. 

An: mean motion difference from computed value. 

■s/as: square root of the semi-major axis of the satellite orbit. 
e s : eccentricity of the satellite orbit. 

T GD , t oc , a /o, a fi , a / 2 - clock correction parameters. 
toe : reference time ephemeris. 

Cus, C uc : amplitude of the sine and cosine harmonic correction term to the 
argument of latitude, respectively. 

C rs , C r( : amplitude of the sine and cosine harmonic correction term to the 
orbit radius, respectively. 

Cis, C ic : amplitude of the sine and cosine harmonic correction term to the 
angle of inclination, respectively. 

U e : longitude of ascending node of orbit plane at weekly epoch. 

Q : rate of the right ascension. 
i: inclination angle at reference time, 
co: argument of perigee, 
idot: rate of inclination angle. 

4.12 SUMMARY 

This chapter takes the satellite position calculated in Chapter 3 and transforms 
it into an earth-centered, earth-fixed coordinate system because this coordinate 
references a fixed position on or above the earth. Since the satellite orbit cannot 
be described perfectly by an elliptic, corrections must be made to the position 
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of the satellite. The information for correction is contained in the ephemeris 
data transmitted by the satellite. This information can be obtained if the GPS 
signal is decoded. The GPS time t at time of transmission needs to be corrected 
for the transit time as well as from the ephemeris data. Obtaining the coarse 
GPS time t c at time of transmission corrected for transit time will be discussed 
in Section 9.7. Finally, the coordinate system of the satellite must be adjusted 
to accommodate the transit time. 
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CHAPTER. FIVE 



GPS C/A Code Signal Structure 



5.1 INTRODUCTION^ 1 ' 2 ) 

In the previous chapters user positions are calculated. In order to perform the 
user position calculation, the positions of the satellites and pseudoranges to 
the satellites must be measured. Many parameters are required to calculate the 
positions of the satellites and they are transmitted in the satellite signals. 

This chapter provides the details associated with the GPS signals. Spilker (L2) 
not only gives a very good discussion on the signal, it also gives the rea- 
sons these signals are selected. The discussion in this chapter is limited to the 
fundamentals of the signals, such that a receiver design can be based on the 
signals. 

There are basically two types of signals: the coarse (or clear)/acquisition 
(C/A) and the precision (P) codes. The actual P code is not directly transmitted 
by the satellite, but it is modified by a Y code, which is often referred to as 
the P(Y) code. The P(Y) code is not available to civilian users and is primarily 
used by the military. In other words, the P(Y) code is classified. The P(Y) code 
has similar properties of the P code. In order to receive the P(Y) code, one must 
have the classified code. Therefore, only the fundamentals of the P code will 
be mentioned in this book. The discussion will be focused on the C/A code. 
In general, in order to acquire the P(Y) code, the C/A code is usually acquired 
first. However, in some applications it is desirable to acquire the P(Y) code 
directly, which is known as direct Y acquisition. 

The radio frequency (RF) of the C/A code will be presented first, then the 
C/A code. The generation of the C/A code and its properties will be presented 
because they are related closely to acquiring and tracking the GPS signals. 
Finally, the data carried by the signals will be presented. The applications of 
the data will be briefly discussed. 
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5.2 TRANSMITTING FREQUENCY! 14 ) 

The GPS signal contains two frequency components: link 1 (LI) and link 2 
(L2). The center frequency of LI is at 1575.42 MHz and L2 is at 1227.6 MHz. 
These frequencies are coherent with a 10.23 MHz clock. These two frequencies 
can be related to the clock frequency as 

LI = 1575.42 MHz = 154 X 10.23 MHz 
L2 = 1227.6 MHz = 120 X 10.23 MHz 

These frequencies are very accurate as their reference is an atomic frequency 
standard. When the clock frequency is generated, it is slightly lower than 10.23 
MHz to take the relativistic effect into consideration. The reference frequency 
is off by® -4.567 x 10 3 Hz, which corresponds to a fraction of -4.4647 x 
10 10 (-4.567 x 10 3 /l 0.23 x 10 6 ). Therefore, the reference frequency used by 
the satellite is 10.229999995433 MHz (10.23 x 10 6 - 4.567 x 10 3 ) rather than 
10.23 MHz. When a GPS receiver receives the signals, they are at the desired 
frequencies. However, the satellite and receiver motions can produce a Doppler 
effect as discussed in Section 3.5. The Doppler frequency shift produced by the 
satellite motion at LI frequency is approximately ±5 KHz. 

The signal structure of the satellite may be modified in the future. However, 
at the present time, the LI frequency contains the C/A and P(Y) signals, while 
the L2 frequency contains only the P(Y) signal. The C/A and P(Y) signals in 
the LI frequency are in quadrant phase of each other and they can be written 
as: 



S L i - A p P(t)D(t) cos(2tt/'i t + <t>) + A c C(t)D(t) sin(2irf 1 t + 0) (5.1) 

where Su is the signal at LI frequency, A p is the amplitude of the P code, P(t) 
- ±1 represents the phase of the P code, D(t) - ±1 represents the data code, 
/ 1 is the LI frequency, <j> is the initial phase, A c is the amplitude of the C/A 
code, C(t ) = ±1 represents the phase of the C/A code. These terms will be 
further discussed in the following sections. In this equation the P code is used 
instead of the P(Y) code. The P(Y), C/A, and the carrier frequencies are all 
phase locked together. 

The minimum power levels of the signals must fulfill the values listed in 
Table 5.1 at the receiver. These power levels are very weak and the spectrum 
is spread, therefore they cannot be directly observed from a spectrum analyzer. 
Even when the signal is amplified to a reasonable power level, the spectrum of 
the C/A code cannot be observed because the noise is stronger than the signal. 

As discussed in Section 3.3, the received power levels at various points on 
the earth are different. The maximum difference is about 2. 1 dB between a point 
just under the satellite and a point tangential to the surface of the earth. In order 
to generate a uniform power over the surface of the earth, the main beam pattern 
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TABLE 5.1 Power Level of GPS Signals 





P 


C/A 


LI 


-133 dBm 


-130 dBm 


L2 


-136 dBm 


-136 dBm* 


‘Presently not ii 


n L2 frequency. 





of the transmitting antenna is slightly weaker at the center to compensate for the 
user at the edge of the beam. The resulting power level versus elevation angle 
is shown in Figure 3.10. The maximum power is -128 dBm, which occurs at 
about 40 degrees. Of course, the receiving antenna pattern also contributes to 
the power level of the receiver. Usually the receiving antenna has a higher gain 
in the zenith direction. This incorporates the ability of attenuating multipath 
but loses gain to signals from lower elevation angles. As discussed in Sections 
3.3 and 3.10, the minimum required beam width of the transmitting antenna 
to cover the earth is 13.87 degrees. The beam width of the antenna® is 21.3 
degrees, which is wider than needed to cover the earth as shown in Figure 5.1. 

If the user is in an aircraft, as long as it is in the main beam of the GPS 
signal and not in the shadow of the earth it can receive the signal. The signals 
generated by the satellite transmitting antenna are right-hand polarized. There- 




FIGURE 5.1 GPS signal main beam. 
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fore, the receiver antenna should be right-hand polarized to achieve maximum 
efficiency. 



5.3 CODE DIVISION-MULTIPLE ACCESS (CDMA) SIGNALS 

A signal S can be written in the following form: 

S = A sin (2-ji ft + 4 >) (5.2) 

where A is the amplitude, / is the frequency, 4> is the initial phase. These three 
parameters can be modulated to carry information. If A is modulated, it is 
referred to as amplitude modulation. Iff is modulated, it is frequency mod- 
ulation. If 4> is modulated, it is phase modulation. 

The GPS signal is a phase- modulated signal with 4> = u,rr; this type of phase 
modulation is referred to as bi-phase shift keying (BPSK). The phase change 
rate is often referred to as the chip rate. The spectrum shape can be described 
by the sine function (sinx/x) with the spectrum width proportional to the chip 
rate. For example, if the chip rate is 1 MHz, the main lobe of the spectrum has 
a null-to-null width of 2 MHz. Therefore, this type of signal is also referred to 
as a spread-spectrum signal. If the modulation code is a digital sequence with a 
frequency higher than the data rate, the system can be called a direct-sequence 
modulated system. 

A code division multiple access (CDMA) signal in general is a spread-spec- 
trum system. All the signals in the system use the same center frequency. The 
signals are modulated by a set of orthogonal (or near-orthogonal) codes. In order 
to acquire an individual signal, the code of that signal must be used to correlate 
with the received signal. The GPS signal is CDMA using direct sequence to 
bi-phase modulate the carrier frequency. Since the CDMA signals all use the 
same carrier frequency, there is a possibility that the signals will interfere with 
one another. This effect will be more prominent when strong and weak signals 
are mixed together. In order to avoid the interference, all the signals should 
have approximately the same power levels at the receiver. Sometimes in the 
acquisition one finds that a cross-correlation peak of a strong signal is stronger 
than the desired peak of a weak signal. Under this condition, the receiver may 
obtain wrong information. 



5.4 P CODE<b2) 

The P code is bi-phase modulated at 10.23 MHz; therefore, the main lobe of the 
spectrum is 20.46 MHz wide from null to null. The chip length is about 97.8 ns 
(1/10.23 MHz). The code is generated from two pseudorandom noise (PRN) 
codes with the same chip rate. One PRN sequence has 15,345,000 chips, which 
has a period of 1.5 seconds, the other one has 15,345,037 chips, and the dif- 
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ference is 37 chips. The two numbers, 15,345,000 and 15,345,037, are relative 
prime, which means there are no common factors between them. Therefore, the 
code length generated by these two codes is 23,017,555.5 (1.5 x 15,345,037) 
seconds, which is slightly longer than 38 weeks. However, the actual length of 
the P code is 1 week as the code is reset every week. This 38-week-long code 
can be divided into 37 different P codes and each satellite can use a differ- 
ent portion of the code. There are a total of 32 satellite identification numbers 
although only 24 of them are in the orbit. Five of the P code signals (33-37) 
are reserved for other uses such as ground transmission. In order to perform 
acquisition on the signal, the time of the week must be known very accurately. 
Usually this time is found from the C/A code signal that will be discussed in 
the next section. The navigation data rate carried by the P code through phase 
modulation is at a 50 Hz rate. 



5.5 C/A CODE AND DATA FORMAT^, 5) 

The C/A code is a bi-phase modulated signal with a chip rate of 1.023 MHz. 
Therefore, the null-to-null bandwidth of the main lobe of the spectrum is 2.046 
MHz. Each chip is about 977.5 ns (1/1.023 MHz) long. The transmitting band- 
width of the GPS satellite in the LI frequency is approximately 20 MHz to 
accommodate the P code signal; therefore, the C/A code transmitted contains 
the main lobe and several sidelobes. The total code period contains 1,023 chips. 
With a chip rate of 1.023 MHz, 1,023 chips last 1 ms; therefore, the C/A code 
is 1 ms long. This code repeats itself every millisecond. The spectrum of a C/A 
code is shown in Figure 5.2. 

In order to find the beginning of a C/A code in the received signal only a 
very limited data record is needed such as 1 ms. If there is no Doppler effect 
on the received signal, then one millimeter of data contains all the 1,023 chips. 
Different C/A codes are used for different satellites. The C/A code belongs to 
the family of Gold codes/ 5 ) which will be discussed in the next section. 

Figure 5.3 shows the GPS data format. The first row shows a C/A code with 
1,023 chips; the total length is 1 ms. The second row shows a navigation data 
bit that has a data rate of 50 Hz; thus, a data bit is 20 ms long and contains 
20 C/A codes. Thirty data bits make a word that is 600 ms long as shown in 
the third row. Ten words make a subframe that is 6 seconds long as shown in 
row four. The fifth row shows a page that is 30 seconds long and contains 5 
subframes. Twenty-five pages make a complete data set that is 12.5 minutes 
long as shown in the sixth row. The 25 pages of data can be referred to as a 
superframe. 

The parameters mentioned in Section 4.10 are contained in the first three 
subframes of a page. If one can receive the information of these three subframes 
from four or more satellites, the user location can be found. Theoretically, one 
can take a minimum of about 18 seconds of data from four satellites and be 
able to calculate the user position. However, the subframes from each satellite 
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FIGURE 5.2 Spectrum of a C/A code. 



will not reach the receiver at the same time. Besides, one does not know when 
the beginning of subframe 1 will be received. A guaranteed way to receive the 
first three subframes is to take 30 seconds (or one page) of data. Thus, one can 
take a minimum of 30 seconds of data and calculate the user position. 



5.6 GENERATION OF C/A CODE(b 2 ' 6 ) 

The GPS C/A signals belong to the family of Pseudorandom noise (PRN) 
codes known as the Gold codes. The signals are generated from the product 
of two 1,023-bit PRN sequence G1 and G2. Both G1 and G2 are generated by 
a maximum-length linear shift register of 10 stages and are driven by a 1.023 
MHz clock. Figure 5.4 shows the G1 and G2 generators. Figure 5.4a shows the 
G1 generator and Figures 5.4b and 5.4c show the G2 generator. Figure 5.4c is 
a simplified notation of Figure 5.4b. 

The basic operating principles of these two generators are similar; therefore, 
only G2 will be discussed in detail. A maximum-length sequence (MLS) gener- 
ator can be made from a shift register with proper feedback. If the shift register 
has n bits, the length of the sequence generated is 2" - 1 . Both shift generators 
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C/A code 1ms 



I C/A code 



i data 20ms 



word 600ms 




in G1 and G2 have 10 bits, thus, the sequence length is 1,023 (2 10 - 1). The 
feedback circuit is accomplished through modulo-2 adders. 

The operating rule of the modulo-2 adder is listed in Table 5.2. When the 
two inputs are the same the output is 0, otherwise it is 1. The positions of the 
feedback circuit determine the output pattern of the sequence. The feedback 
of G1 is from bits 3 and 10 as shown in Figure 5.4a and the corresponding 
polynomial can be written as Gl: 1 + x 3 + x'°. The feedback of G2 is from 
bits 2, 3, 6, 8, 9, 10 as shown in Figure 5.4b and the corresponding polynomial 
is G2: 1 + x 2 + x 3 + x 6 + x 8 + x 9 + x 10 . 

In general, the output from the last bit of the shift register is the output of 
the sequence as shown in Figure 5.4a. Let us refer to this output as the MLS 
output. However, the G2 generator does not use the MLS output as the output. 
The output is generated from two bits which are referred to as the code phase 
selections through another modulo-2 adder as shown in Figures 5.4b and c. This 
G2 output is a delayed version of the MLS output. The delay time is determined 
by the positions of the two output points selected. 
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(c) Simplified notation of G2. 

FIGURE 5.4 Gl, G2 maximum-length sequence generators. 



Figure 5.5 shows the C/A code generator. Another modulo-2 adder is used 
to generate the C/A code, which uses the outputs from Gl and G2 as inputs. 
The initial values of the two shift registers Gl and G2 are all l’s and they must 
be loaded in the registers first. The satellite identification is determined by the 



TABLE 5.2 Modulo-2 Addition 



Input 1 


Input 2 


Output 


0 


0 


0 


0 


1 


1 


1 


0 


1 


1 


1 


0 
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G1 generator 




FIGURE 5.5 C/A code generator. 



two output positions of the G2 generator. There are 37 unique output positions. 
Among these 37 outputs, 32 are utilized for the C/A codes of 32 satellites, 
but only 24 satellites are in orbit. The other five outputs are reserved for other 
applications such as ground transmission. 

Table 5.3 lists the code phase assignments. In this table there are five 
columns and the first column gives the satellite ID number, which is from 1 
to 32. 

The second column gives the PRN signal number; and it is from 1 to 37. 
It should be noted that the C/A codes of PRN signal numbers 34 and 37 are 
the same. The third column provides the code phase selections that are used to 
form the output of the G2 generator. The fourth column provides the code delay 
measured in chips. This delay is the difference between the MLS output and the 
G2 output. This is redundant information of column 3, because once the code 
phase selections are chosen this delay is determined. The last column provides 
the first 10 bits of the C/A code generated for each satellite. These values can be 
used to check whether the generated code is wrong. This number is in an octal 
format. 

The following example will illustrate the use of the information listed in 
Table 5.3. For example, in order to generate the C/A code of satellite 19, the 
3 and 6 tabs must be selected for the G2 generator. With this selection, the G2 
output sequence is delayed 471 chips from the MLS output. The last column is 
1633, which means 1 110011 Oil in binary form. If the first 10 bits generated 
for satellite 19 do not match this number, the code is incorrect. 
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TABLE 5.3 Code Phase Assignments 



Satellite ID 


GPS PRN 


Code Phase 


Code Delay 


First 10 Chips 


Number 


Signal Number 


Selection 


Chips 


C/A Octal 


1 


1 


2 


© 


6 


5 


1440 


2 


2 


3 


© 


7 


6 


1620 


3 


3 


4 


© 


8 


7 


1710 


4 


4 


5 


© 


9 


8 


1744 


5 


5 


1 


© 


9 


17 


1133 


6 


6 


2 


© 


10 


18 


1455 


7 


7 


1 


© 


8 


139 


1131 


8 


8 


2 


© 


9 


140 


1454 


9 


9 


3 


© 


10 


141 


1626 


10 


10 


2 


© 


3 


251 


1504 


11 


11 


3 


© 


4 


252 


1642 


12 


12 


5 


© 


6 


254 


1750 


13 


13 


6 


© 


7 


255 


1764 


14 


14 


7 


© 


8 


256 


1772 


15 


15 


8 


© 


9 


257 


1775 


16 


16 


9 


© 


10 


258 


1776 


17 


17 


1 


© 


4 


469 


1156 


18 


18 


2 


© 


5 


470 


1467 


19 


19 


3 


© 


6 


471 


1633 


20 


20 


4 


© 


7 


472 


1715 


21 


21 


5 


© 


8 


473 


1746 


22 


22 


6 


© 


9 


474 


1763 


23 


23 


1 


© 


3 


509 


1063 


24 


24 


4 


© 


6 


512 


1706 


25 


25 


5 


© 


7 


513 


1743 


26 


26 


6 


© 


8 


514 


1761 


27 


27 


7 


© 


9 


515 


1770 


28 


28 


8 


© 


10 


516 


1774 


29 


29 


1 


© 


6 


859 


1127 


30 


30 


2 


© 


7 


860 


1453 


31 


31 


3 


© 


8 


861 


1625 


32 


32 


4 


© 


9 


862 


1712 




33 


5 


© 


10 


863 


1745 




34* 


4 


© 


10 


950 


1713 




35 


1 


© 


7 


947 


1134 




36 


2 


© 


8 


948 


1456 




37* 


4 


© 


10 


950 


1713 



*34 and 37 have the same C/A code. 

**GPS satellites do not transmit these codes; they are reserved for other uses. 
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A computer program (p5_l) is listed at the end of this chapter to generate 
both the MLS and G2 output sequences. The program takes columns 3 and 4 
of Table 5.3 as inputs and checks the time delay. If the correct data are used 
as inputs, the output will show “OK,” otherwise, it will show “not match.” 

A program (p5_2) can be used to generate the C/A code. The program is an 
extension of the program (p5_l) to include the two maximum-length sequence 
generators. In the program, the delay time listed in Table 5.3 is used as input 
to generate the G2 signal rather than using the code phase selections in column 
3. The first 10 bits of the generated C/A code should be compared with the 
result listed in the last column of Table 5.3. 



5.7 CORRELATION PROPERTIES OF C/A CODE(L2) 

One of the most important properties of the C/A codes is their correlation result. 
High autocorrelation peak and low cross-correlation peaks can provide a wide 
dynamic range for signal acquisition. In order to detect a weak signal in the 
presence of strong signals, the autocorrelation peak of the weak signal must be 
stronger than the cross-correlation peaks from the strong signals. If the codes 
are orthogonal, the cross correlations will be zero. However, the Gold codes 
are not orthogonal but near orthogonal, implying that the cross correlations are 
not zero but have small values. 

The cross correlation of the Gold code is listed in Table 5.4.0) 



TABLE 5.4 


Cross Correlation of Gold Code 






Code Period 


Number of Shift Normalized Cross 

Register Stages Correlation Level 


Probability of 
Level 






r 2(" + 1 )/2 + l 

, P 


0.25 


P = 2" - 1 


n = odd \ 


1 1 
1 P 


0.5 






2 {n+2)/2_ | 


0.24 






L P 




1 


2 (n+2)/2 + i 

I P 


0.125 


P = 2" - 1 


n = even < 


1 1 
1 P 


0.75 




1 


| 2 (» + 2)/2_ 1 


0.125 
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(b) Cross correlation of satellites 19 and 31. 
FIGURE 5.6 Auto and cross correlation of C/A code. 



For the C/A code n - even = 10, thus, P - 1023. Using the relations in the above 
table, the cross-correlation values are: -65/1023 (occurrence 12.5%), -1/1023 
(75%), and 63/1023 (12.5%). The autocorrelation of the C/A codes of satellite 
19 and the cross correlation of satellites 19 and 31 are shown in Figures 5.6a 
and 5.6b respectively. These satellites are arbitrarily chosen. 

In Figure 5.6a, the maximum of the autocorrelation peak is 1023, which 
equals the C/A code length. The position of the maximum peak is deliberately 
shifted to the center of the figure for a clear view. The rest of the correlation 
has three values 63, -1, and -65. The cross-correlation shown in Figure 5.6b 
also has three values 63, -1, -65. 

These are the values calculated by using equations in Table 5.4. The difference 
between the maximum of the autocorrelation to the cross correlation determines 
the processing gain of the signal. In order to generate these figures, the outputs 
from the C/A code generator must be 1 and - 1, rather than 1 and 0. The mathe- 
matical operation to generate these figures will be discussed in the Section 7.7. 



5.8 NAVIGATION DATA BITS* 2 ' 3 ' 7 ) 

The C/A code is a bi-phase coded signal which changes the carrier phase 
between 0 and -k at a rate of 1.023 MHz. The navigation data bit is also bi- 





